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Abstract 


This  study  proved  that  transpiration  cooling  provides  a  better  cooling  scheme  than 
regenerative  cooling  for  long  operating  duration  liquid-ftieled  rocket  engine  nozzles.  This 
proof  was  made  on  the  basis  of  maximum  wall  temperature.  This  study  compared 
transpiration  cooling  to  regenerative  cooling  in  the  throat  region  of  the  Space  Shuttle 
Main  Engine  Main  Combustion  Chamber.  The  study  also  analyzed  the  effects  of  porosity, 
solid  thermal  conductivity,  and  porous  sphere  size  on  a  porous  wall  made  of  packed 
spheres.  The  transpiration  cooled  nozzle  operated  35%  cooler  than  a  regeneratively 
cooled  nozzle,  but  the  temperature  gradient  at  the  hot  gas  surface  was  72  times  greater 
than  the  regeneratively  cooled  nozzle. 
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NUMERICAL  STUDY  OF  A  TRANSPIRATION  COOLED  ROCKET  NOZZLE 


I.  Introduction 


1.1  Background 

Future  space  systems  are  looking  toward  higher  thrust,  longer  operating  duration, 
more  reliable,  more  reusable  systems.  To  increase  thrust  efficiency,  combustion  chambers 
must  operate  at  higher  chamber  pressures  and  temperatures.  This  requires  sufficient 
cooling  to  maintain  the  combustion  chamber  material  within  acceptable  temperature  limits 
either  to  prevent  failure  or  increase  cycle  life.  The  peak  heat  flux  in  a  rocket  nozzle  occurs 
near  the  throat,  therefore  the  true  measure  of  a  cooling  system  is  how  well  it  operates  in 
the  throat  region  (Hill  and  Peterson,  1992:550). 

The  currently  accepted  norm  for  cooling  the  nozzles  of  long  operating  duration, 
liquid-fueled  rocket  engines  is  regenerative  cooling  (Mueggenburg  and  others,  1992:2). 
Regenerative  cooling,  or  forced  convection  cooling,  passes  cold  fuel  at  high  velocities 
through  small  channels  imbedded  in  the  nozzle  wall  to  absorb  the  heat  from  the  hot  wall 
and  to  preheat  the  fuel  (Hill  and  Peterson,  1992:541).  Another  method  to  cool  the  rocket 
nozzle  is  transpiration  cooling.  Transpiration  cooling  is  the  injection  of  a  cooling  fluid 
through  a  porous  wall  material  over  a  relatively  large  surface  area. 
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1.2  Thesis  Statement 


Transpiration  cooling  provides  a  better  cooling  scheme  than  regenerative  cooling 
for  a  liquid-fueled,  long  operating  duration  rocket  nozzle. 

1.3  Justification 

The  most  technologically  advanced  example  of  regenerative  cooling  used  in  a 
reusable,  liquid-fueled,  high  chamber  pressure  rocket  engine  is  the  Space  Shuttle  Main 
Engine  (SSME).  Even  though  regenerative  cooling  was  originally  chosen  for  the  SSME, 
after  several  operational  flights  it  was  discovered  that  regenerative  cooling  does  not 
provide  enough  cooling  for  the  SSME  to  meet  cycle  life  requirements.  Thermally  induced 
stresses  caused  surface  defects  which  increased  the  number  of  between-launch  service 
repairs.  Two  identified  failure  modes  resulted  in  cracks  in  the  surface  of  the  SSME  main 
combustion  chamber  (MCC)  (Murphy  and  others,  1986:2;  Quentmeyer,  1990:1). 

Regenerative  cooling  also  creates  other  problems.  One  problem  is  the  pressure 
drop  that  occurs  in  the  small  coolant  channels.  Regenerative  cooling  requires  high 
velocity  coolant  flow  through  very  small  channels  to  increase  heat  transfer  enough  to  cool 
the  hot  wall.  This  causes  a  large  pressure  drop  through  the  coolant  channels.  This 
pressure  drop  creates  one  of  two  problems.  Either  the  fuel  turbopump  requires  a  higher 
discharge  pressure  to  maintain  the  same  combustion  chamber  pressure,  or  the  chamber 
pressure  must  decrease  to  make  up  for  the  pressure  drop.  Furthermore,  higher  turbopump 
discharge  pressures  cause  a  decrease  in  life  of  the  turbopump  (May  and  others,  1990: 1). 

Another  problem  with  regenerative  cooling  is  that  it  requires  very  thin  hot  gas  wall 
thicknesses.  This  is  done  to  prevent  steep  temperature  gradients  that  cause  thermal  strain 
and  decreased  cycle  life.  Unfortunately,  the  desired  wall  thicknesses  approach  the 
tolerances  of  conventional  machining  processes  (Mueggenburg  and  others,  1992:2). 
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Because  of  these  problems,  the  National  Aeronautics  and  Space  Administration 
(NASA)  and  the  Air  Force  Advanced  Launch  System  (ALS)  began  looking  for  high 
reliability,  low  cost  cooling  schemes  for  proposed  future  booster  designs.  A  study  in  1990 
summarized  several  new  cooling  techniques  and  pointed  out  transpiration  cooling  as  a 
promising  new  method  (Quentmeyer,  1990:5-6). 

Work  has  been  done  at  the  Air  Force  Institute  of  Technology  (AFIT)  to 
experimentally  explore  the  effects  of  transpiration  cooling  on  heat  transfer  to  the  wall  of  a 
supersonic  nozzle  (Lenertz,  1994;  Keener,  1994;  Chen,  1995).  The  most  important  result 
of  these  studies  was  a  simple  mathematical  expression  that  quantifies  the  decrease  in  heat 
transfer  to  the  hot  gas  side  wall  due  to  the  injection  of  the  coolant  (Chen,  1995:1). 

1.4  Purpose 

The  purpose  of  this  thesis  is  to  prove  that  transpiration  cooling  provides  a  better 
cooling  scheme  than  regenerative  cooling  for  a  liquid-fueled,  long  operating  duration 
rocket  nozzle.  To  prove  this  fact,  the  author  chose  maximum  wall  temperature  and 
maximum  temperature  gradient  as  the  figures  of  merit  for  cooling  performance.  A  cooling 
scheme  that  provides  the  lowest  wall  temperature  with  the  lowest  temperature  gradient 
would  be  considered  the  best  cooling  scheme. 

This  study  also  characterizes  the  relative  importance  of  different  geometric  and 
thermodynamic  properties  of  the  wall  material.  These  properties  are  the  porosity,  the  void 
fraction  of  the  porous  material,  and  the  thermal  conductivity,  the  material's  ability  to 
transfer  heat  through  the  solid.  Another  parameter  of  interest  is  the  blowing  ratio.  The 
blowing  ratio  is  defined  as  the  mass  flow  per  unit  area  of  injected  fluid  divided  by  the  mass 
flow  per  unit  area  of  free  stream  hot  gas.  Or  in  mathematical  terms,  the  blowing  ratio  is 
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BR  = 


(1) 


Pi^i 
p  u 

g 

where  BR  is  the  blowing  ratio,  p,  is  the  density  of  the  injected  coolant,  m,  is  the  uniform 
velocity  of  the  injected  coolant  at  the  porous  surface,  p^  is  the  density  of  the  free  stream 
hot  gases,  and  is  the  uniform  free  stream  velocity  of  the  hot  gases  (Lenertz,  1994:2). 

The  pressure  difference  across  the  porous  nozzle  wall  of  a  transpiration  cooled  nozzle 
determines  the  blowing  ratio.  By  studying  the  relative  importance  of  varying  these  three 
parameters,  the  design  engineer  will  be  able  to  select  the  appropriate  porosity,  material  to 
construct  the  porous  nozzle  wall,  and  blowing  ratio  for  a  given  problem. 

1.5  Methodology 

The  author  developed  a  finite  difference  numerical  model  of  the  heat  transfer 
through  the  porous  wall.  This  model  treats  the  transpiration  coolant  flow  and  heat 
transfer  as  one  dimensional.  This  assumption  is  justified  through  the  advancement  of 
platelet  technology.  Platelet  technology  involves  photo-etching  channels  into  very  thin 
layers  of  material  and  bonding  these  layers  together  to  form  a  composite  structure.  These 
tiny  chaimels  provide  accurate  metering  of  the  flow  to  meet  either  a  specific  pressure  or 
heat  flux  profile  without  two  or  three  dimensional  coolant  flow  effects  (Mueggenburg  and 
others,  1992;  1-2).  The  one  dimensional  heat  transfer  assumption  is  justified  due  to  the 
relative  magnitudes  of  the  temperature  gradients.  The  temperature  gradient  along  the 
nozzle  wall  does  not  approach  the  4x10^  K/m  gradient  that  occurs  within  the  nozzle  wall, 
therefore  two  dimensional  effects  can  be  neglected. 

The  author  chose  to  represent  the  porous  nozzle  wall  as  a  packed  bed  of  spheres 
with  the  transpiration  coolant  constrained  to  flow  in  only  one  direction.  Since  the 
transpiration  coolant  follows  a  tortuous  path  through  both  the  platelet  wall  and  the  packed 
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bed  of  spheres  wall  and  the  coolant  flow  is  constrained  to  flow  in  only  one  direction,  the 
packed  bed  of  spheres  model  is  valid.  Five  different  regular  geometrical  arrangements  of 
spherical  particles  were  used  to  represent  the  porous  material  construction.  This  allowed 
computation  of  porosity  and  solid  surface  area  exposed  to  fluid  flow.  The  solid  surface 
area  is  critical  in  linking  the  heat  transfer  between  the  solid  and  the  coolant  as  it  flows 
through  the  porous  solid.  In  addition,  the  use  of  spheres  allowed  the  author  to  use  the 
heat  transfer  equations  for  heat  transfer  from  a  gas  to  a  packed  bed  of  spheres. 

Since  the  SSME  MCC  is  the  most  technologically  advanced,  operational  example 
of  regenerative  cooling,  the  author  chose  it  as  the  standard  against  which  all  transpiration 
cases  should  be  judged.  The  SSME  MCC  nozzle  and  coolant  channel  dimensions  were 
used  in  the  heat  transfer  equations.  Because  hydrogen  is  already  used  as  the  cooling  fluid 
in  the  SSME  MCC  at  the  throat  region,  hydrogen  was  chosen  as  the  transpiration  coolant 
fluid.  Additionally,  hydrogen  would  be  a  good  coolant  for  transpiration  cooling  because  it 
has  a  low  molecular  weight,  which  is  desirable  for  a  higher  engine  specific  impulse 
(Bowman  and  others,  1994:20).  At  nozzle  pressures  and  temperatures,  the  hydrogen 
coolant  is  a  supercritical  fluid  that  behaves  like  an  ideal  gas.  The  hydrogen  is  neither  a 
liquid  or  a  gas,  it  is  in  a  state  between  the  two.  Throughout  this  thesis,  the  terms  fluid, 
liquid  and  gas  are  used  interchangeably  to  refer  to  the  transpiration  coolant. 

Five  different  thermal  conductivities  representing  five  different  solid  materials  were 
used  in  the  parameter  study.  The  values  of  these  thermal  conductivities  represent  the 
range  from  insulators  to  very  good  thermal  conductors.  Five  different  pressure  differences 
were  used  to  obtain  a  spread  of  blowing  ratios.  An  additional  parameter  in  the  porous 
material  is  the  sphere  size  used  to  make  the  porous  solid.  There  were  three  sphere  radius 
sizes  chosen. 
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The  numerical  model  was  built  in  a  modular  design.  Each  subroutine  was  tested 
independently.  The  solid  and  liquid  subroutines  were  validated  against  known  exact 
solutions.  Once  these  subroutines  were  validated,  they  were  combined  to  create  the 
complete  model.  The  model  was  run  comparing  each  of  the  porosities  against  each  of  the 
thermal  conductivities  at  various  pressure  differences  to  generate  a  spread  of  blowing 
ratios.  Finally,  the  maximum  temperatures  and  maximum  temperature  gradients  were 
plotted  against  blowing  ratio  to  observe  trends  as  porosity  and  thermal  conductivity  vary. 
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II.  Numerical  Model 


This  chapter  provides  descriptions  of  the  equations  used  in  the  numerical  model. 
The  chapter  begins  with  an  overview  of  the  heat  transfer  processes  within  a  porous  wall 
and  a  statement  of  the  energy  equations  for  the  porous  solid  and  liquid  transpiration 
coolant.  The  next  two  sections  describe  the  equations  that  govern  the  boundary  heat  flux 
conditions.  The  fourth  section  in  this  chapter  describes  the  heat  transfer  that  takes  place 
between  the  porous  solid  and  liquid  coolant.  The  fifth  section  describes  the  source  of  the 
properties  used  in  the  model.  The  final  section  of  this  chapter  gives  a  brief  overview  of 
the  program  structure.  The  model  was  programmed  in  FORTRAN  and  the  code  is 
provided  in  Appendix  B. 

2. 1  Transpiration  Cooling  Overview 

The  various  heat  and  energy  transfer  mechanisms  involved  in  transpiration  cooling 
are  illustrated  in  Figure  1 .  The  two  mechanisms  of  heat  transfer  used  in  this  study  are 
conduction  and  convection.  One  dimensional  conduction  heat  transfer  is  described  by 
Fourier's  Law 


q 


n _ 


dT 
k  — 
dx 


(2) 


where  q"  is  the  heat  flux  due  to  conduction  in  watts  per  square  meter,  k  is  the  thermal 

dT 

conductivity  of  the  material,  and  —  is  the  temperature  gradient  within  the  material 

dx 

(Incropera  and  DeWitt,  1990:4).  Convection  heat  transfer  between  a  solid  of  temperature, 
7^,  and  a  liquid  of  temperature,  is  described  by  Newton's  Law  of  Cooling 
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g"  =  h{T,-T,) 


(3) 


where  q"  is  the  heat  flux  due  to  convection  in  watts  per  square  meter,  and  h  is  the  heat 
transfer  coefficient  in  watts  per  square  meter  per  degree  Kelvin  (Incropera  and  DeWitt, 
1990:8).  Energy  is  also  transferred  through  the  transpiration  coolant  by  the  work  done  by 
the  pressure  forces  on  the  fluid  and  the  advection  of  kinetic  and  internal  energy  by  the 
fluid  (Incropera  and  DeWitt,  1990:326-334). 

The  heat  transfer  and  energy  mechanisms  are  combined  to  form  the  energy 
equations  for  the  solid  and  the  liquid.  Appendix  A  contains  the  derivations  of  the  energy 
equations,  but  the  final  forms  are  provided  here  for  the  reader.  The  solid  energy  equation 
assumes  (1)  the  only  mode  of  heat  transfer  through  the  solid  is  by  conduction,  and  (2)  the 
only  loss  of  energy  is  to  the  transpiration  coolant.  The  energy  equation  for  the  porous 
solid  used  in  this  model  is  given  in  Equation  (4). 


,.d^T,  KA-(T.-T,)  ^ 

dx^  Actx  '  ’’’  dt 


(4) 


The  first  term  in  the  solid  energy  equation  represents  the  conduction  through  the  solid. 
The  second  term  describes  the  heat  transfer  from  the  solid  to  the  liquid,  where  A*  is  the 
surface  area  of  the  solid  exposed  to  the  transpiration  coolant  and  A  is  the  cross  sectional 
area  of  the  nozzle  wall  exposed  to  the  transpiration  coolant  flow.  The  term  on  the  right 
hand  side  of  the  energy  equation  is  the  transient  term  and  is  only  required  for  obtaining 
convergence  to  the  steady  state  solution. 
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The  words  above  the  equations  indicate  the  direction  of  heat  transfer,  e.g.  Solid-Coolant  indicates  heat  transfer 
from  the  solid  to  the  transpiration  coolant 


Figure  1.  Heat  Transfer  Mechanisms  Within  a  Porous  Solid 

The  liquid  energy  equation  assumed;  (1)  the  only  work  done  on  the  fluid  by 
external  forces  is  by  pressure  forces,  (2)  the  only  means  of  energy  transfer  across  the 
control  volume  is  by  conduction  and  advection,  (3)  the  temperature  gradient  term 
dominates  the  heat  transfer  and  the  kinetic  energy  term  can  be  neglected,  (4)  the  coolant 
behaves  like  an  ideal  gas,  and  (5)  the  heat  transfer  from  the  porous  solid  to  the  coolant  is 


9 


the  only  energy  gain.  The  energy  equation  for  the  liquid  transpiration  coolant  used  in  this 
model  is  given  in  Equation  (5). 


‘  dx^  A,dx 


-p{Uf 


dT 


I  _ 


Pi 


dx 


=  0 


(5) 


The  first  term  in  the  liquid  energy  equation  is  the  conduction  through  the  liquid.  The 
second  term  describes  the  amount  of  energy  that  is  transferred  fi'om  the  solid  to  the  liquid 
and  is  considered  an  energy  gain.  The  third  term  describes  the  advection  of  energy 
through  the  liquid  and  is  considered  an  energy  loss. 

The  differential  energy  equations  were  converted  to  explicit  finite  difference  forms 
so  they  could  be  solved  numerically.  The  TSOLED  subroutine  solves  for  the  solid 
temperature  profile  and  the  TLIQUID  subroutine  calculates  the  liquid  temperature  profile. 
Before  the  energy  equations  could  be  solved  numerically,  the  boundary  conditions  on  the 
hot  gas  and  cold  gas  sides  must  be  established. 

2.2  Hot  Gas  Side  Heat  Transfer 

The  heat  transfer  on  the  hot  combustion  gas  side  of  the  nozzle  wall  comes  from 
convection  and  radiation.  Because  the  heat  transfer  correlation  used  includes  radiation 
effects,  the  author  did  not  add  the  radiation  boundary  condition  to  the  energy  equation 
(Bartz,  1965:16-17).  Newton's  law  of  cooling  states  that  the  heat  flux  due  to  convection 
from  a  hot  gas  to  a  solid  surface  takes  the  form 

q"  =  h^{T^-T^)  (6) 
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where  is  the  heat  transfer  coefficient,  7^  is  the  fluid  temperature,  and  is  the  solid 

wall  surface  temperature  at  the  point  of  interest.  For  high  speed  flow,  the  fluid 
temperature,  7^,  is  replaced  with  the  adiabatic  wall  temperature,  due  to  viscous 

heating  in  the  turbulent  boundary  layer.  This  adiabatic  wall  temperature  can  be  found 
fi-om  the  recovery  factor,  R ,  represented  by 

y’  _  T* 

7?  =  Pr^  =  — L  (7) 

where  7^^  is  the  stagnation  temperature  of  the  fluid,  and  7^  is  the  free  stream  temperature 

of  the  fluid  (Kreith  and  Bohn,  1986:586).  The  Prandtl  number,  Pr,  of  the  fluid  is 
evaluated  at  the  film  temperature,  7). 


7}  =  +  0. 23(T,  -  t;;^  ) + 0. 1 9(  t;,,  -  )  (8) 


(Hill  and  Peterson,  1992:550).  The  FINDTAW  subroutine  solves  for  the  adiabatic  wall 
temperature. 

The  last  term  needed  is  the  heat  transfer  coefficient,  .  An  expression  derived  by 

Bartz  for  the  heat  transfer  coefficient  fi'om  the  hot  combustion  gases  to  the  nozzle  wall 
without  transpiration  cooling  effects  is  given  in  Equation  (9). 


0.026 
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where  represents  heat  transfer  without  blowing,  Z),  and  are  the  nozzle  throat 

diameter  and  area  respectively,  is  the  radius  of  curvature  of  the  nozzle  in  the  plane 
that  contains  the  nozzle  axis,  is  the  mass  flow  rate  of  the  combustion  products,  4  is 

the  nozzle  cross  sectional  area  of  interest,  and  the  0  subscript  implies  that  these  properties 
are  evaluated  at  stagnation  conditions  (Bartz,  1965:35).  The  cr  term  is  given  by 


- 

- 

0.8-02a) 

1  T 
'■ 

fl+^“^M0+  ^ 
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_ 1 

2  T 
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^2  j  2_ 
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where  y  is  the  ratio  of  specific  heats  of  the  combustion  products,  M  is  the  local  Mach 
number  of  the  combustion  gases,  and  a  is  the  exponent  of  the  viscosity-temperature 
relation,  p  oc  T®.  For  diatomic  gases,  <0  =  0.6  (Hill  and  Peterson,  1992:550).  Since  the 
combustion  products  are  not  all  diatomic  gases,  the  author  plotted  the  viscosity  of  the 
combustion  products  at  the  nozzle  throat  conditions  against  temperature  and  obtained  the 
value  of  <y  =  0.875.  The  approximate  error  in  using  this  correlation  is  2%  but  can  be  as 
high  as  30%  near  the  throat  (Bartz,  1965:59). 

Lenertz  and  Chen  (1994;  1995)  discovered  that  there  is  a  decrease  in  heat  flux 
from  a  hot  gas  to  the  nozzle  wall  if  transpiration  cooling  is  used  through  the  nozzle  wall. 
An  objective  in  transpiration  cooling  is  to  maximize  the  heat  reduction,  while  minimizing 
the  disturbance  to  the  primary  flow  boundary  layer  (Sucec,  1985:821).  If  the  blowing  is 
too  great,  the  boundary  layer  grows  and  the  effective  nozzle  throat  and  nozzle  exit  areas 
are  decreased  causing  a  decrease  in  performance.  Even  at  low  blowing  ratios.  Keener 
(1994:65)  found  there  is  a  47%  increase  in  the  boundary  layer  with  transpiration  cooling 
although  much  of  this  increase  was  due  to  the  rough  surface.  The  range  of  blowing  ratios 
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considered  by  Chen  was  from  -0.0002  to  0.01 17,  where  negative  values  indicate  suction 
and  positive  values  indicate  blowing.  Lenertz  and  Chen  (1994;  1995)  used  a  two 
dimensional  Mach  2  nozzle  in  the  AFIT  shock  tube  to  measure  the  effect  of  transpiration 
cooling  on  heat  transfer.  Chen  created  a  mathematical  expression  that  related  the  blowing 
heat  transfer  coefficient  to  the  non-blowing  heat  transfer  coefficient  in  Equation  (11). 


-^  =  (\-3S.0*BR) 

^go 


(11) 


where  is  the  heat  transfer  including  blowing  effects,  is  heat  transfer  without 

blowing,  and  BR  is  the  blowing  ratio  (Chen,  1995:4. 1 1).  Since  the  maximum  blowing 
ratio  valid  for  Chen's  correlation  was  0.01 17  (1995:1),  the  author  chose  to  limit  the  scope 
of  blowing  ratios  to  less  than  0.01.  The  HOTGAS  subroutine  calculates  the  hot  gas  side 
heat  transfer  coefficient  including  transpiration  effects. 

2.3  Cold  Gas  Side  Heat  Transfer 

Like  the  hot  gas  side,  the  heat  transfer  on  the  cold  gas  side  is  defined  by 
convection.  The  heat  flux  from  convection  from  a  hot  solid  to  a  cold  gas  takes  the  form 

q"  =  K{T^c-T,)  (12) 


where  is  the  heat  transfer  coefficient,  7^  is  the  bulk  fluid  temperature,  and  7^^  is  the 
solid  wall  surface  temperature.  Since  7^^  is  not  known  initially,  the  author  chose  to  use 
an  expression  for  the  heat  transfer  coefficient  that  accounts  for  large  temperature 
differences  between  7^  and  7^^ ,  and  the  cold  fluid  is  in  turbulent  conditions.  The 
expression  used  is  given  in  Equation  (13)  (Incropera  and  DeWitt,  1990:496). 
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(13) 
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\0.14 
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Re/5 
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where  is  the  coolant  channel  hydraulic  diameter,  Re^^  is  the  Reynolds  number  defined 
as 


Re 


(14) 


where  is  the  mass  flow  rate  of  the  coolant  in  the  coolant  channel,  and  Per^  is  the 

coolant  channel  wetted  perimeter.  In  Equation  (13)  all  properties  are  evaluated  at  the 
coolant  bulk  temperature,  7^ ,  except  which  is  evaluated  at  the  solid  surface 

temperature,  7^^ .  The  approximate  error  in  using  this  correlation  is  25%  (Incropera  and 

DeWitt,  1990:496).  The  HCOOLSIDE  subroutine  calculates  the  heat  transfer  coefficient, 

K- 

2.4  Packed  Bed  of  Spheres  Heat  Transfer 

The  only  remaining  unknown  in  the  energy  equations  is  the  heat  transfer  coefficient 
for  convection  between  the  porous  solid  and  the  transpiration  coolant.  Since  the  porous 
solid  is  assumed  to  be  made  of  a  packed  bed  of  very  small  spheres,  an  expression  for  the 
heat  transfer  coefficient,  h^^,  is  given  in  Equation  (15)  (Whitaker,  1972:366-368). 

//  (o.5Re^  +  0.2Re,^)pr^  (15) 

2r  s  ^  ^ ' 
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where  is  the  thermal  conductivity  of  the  coolant,  is  the  sphere  radius  of  the  porous 
solid  material,  Pr  is  the  Prandtl  number  of  the  fluid,  and  Re,,^  is  the  Reynolds  number 

defined  by 


2rrh, 

Re,  = - 


(16) 


where  ih^  is  the  mass  flow  rate  of  the  coolant  through  the  porous  solid,  and  s  is  the 
porosity  of  the  solid.  All  liquid  properties  in  Equation  (15)  are  evaluated  at  the  average 
between  the  bulk  liquid  temperatures  entering  and  exiting  the  porous  wall  (Whitaker, 
1972:366-368). 

Since  the  flow  velocity  of  the  coolant  is  relatively  slow  in  the  porous  solid,  the 
coolant  is  assumed  to  be  incompressible.  Darcy's  law  relates  the  mass  flow  rate  of  an 
incompressible  liquid  to  the  pressure  gradient  in  the  liquid  in  Equation  (17)  (Brennan  and 
Kroliczek,  1979:25). 


dP  _ 

dx  Kp,A, 


(17) 


In  Equation  (17),  —  is  the  pressure  gradient  through  the  porous  solid,  rh;  is  the  mass 
dx 

flow  rate  of  the  coolant  through  the  porous  solid,  /J■^  and  p,  are  the  viscosity  and  density 
of  the  coolant,  and  K  is  the  permeability  of  the  solid.  For  a  packed  bed  of  spheres,  the 
permeability,  K,  is  defined  in  Equation  (18)  (Brennan  and  Kroliczek,  1979: 118). 


K  = 


150(1- ff)^ 


(18) 
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Since  the  heat  transfer  coefficient,  h^i ,  is  evaluated  at  the  average  temperature,  the 
viscosity  and  density  should  be  evaluated  at  the  same  average  temperature.  Therefore  all 
the  terms  are  constant  throughout  the  porous  solid  and  Equation  (17)  simplifies  to 

L  Kp, 

where  is  the  coolant  line  pressure,  P,  is  the  nozzle  static  throat  pressure,  and  L  is  the 
thickness  of  the  porous  wall.  This  equation  can  be  rearranged  to  obtain  /w,  in  Equation 
(20). 


m 


_{P,~P,)Kp, 

LPi 


(20) 


The  HCOOL  subroutine  calculates  the  heat  transfer  coefficient,  h^, . 

2.5  Properties 

To  calculate  the  heat  transfer  at  the  hot  gas  side  of  the  porous  wall,  the  author 
needed  thermodynamic  properties  for  the  hot  combustion  gases  at  various  temperatures. 
The  combustion  product  properties  came  from  Mr.  Dan  Griesen  at  AEROJET  corporation 
in  Sacramento,  California.  The  properties  were  generated  by  a  computer  code  from 
NASA  called  TRAN72.  Mr.  Griesen  provided  the  author  with  the  output  from  a  6: 1  fuel 
to  oxidizer  combustion  mixture  at  20.68  MPa  and  3656.7  K  combustion  temperature. 
These  conditions  represent  100%  thrust  of  a  SSME.  Curve  fits  were  made  for  these  hot 
gas  properties.  The  HOTGASPROP  subroutine  contains  the  curve  fits  for  the  combustion 
product  properties  (Griesen,  1995). 
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Likewise,  the  author  needed  various  thermodynamic  properties  for  hydrogen  at 
varying  pressures  and  temperatures  to  calculate  the  packed  bed  of  sphere  heat  transfer  and 
the  cold  gas  side  heat  transfer.  The  hydrogen  coolant  properties  came  from  Mr.  Don 
Rousar  at  AEROJET  corporation  in  Sacramento,  California.  They  are  from  an  internal 
report  generated  by  AEROJET  for  use  in  their  calculations.  Curve  fits  were  made  for  the 
pertinent  hydrogen  properties.  The  H2PROP  subroutine  contains  the  curve  fits  for  the 
hydrogen  properties  (Rousar,  1995). 

The  SSME  MCC  geometrical  dimensions  and  the  NARloy-Z  properties  were 
obtained  from  Ms.  Kim  Parks  at  Marshall  Space  Flight  Center  in  Huntsville,  Alabama. 

She  provided  technical  drawings  of  the  SSME  MCC  and  plots  of  NARloy-Z  thermal 
conductivity.  She  also  provided  a  computer  generated  output  of  several  SSME  key 
properties  such  as,  coolant  mass  flow  rate  through  the  coolant  channels,  and  combustion 
product  mass  flow  rate  through  the  nozzle  throat  (Parks,  1995). 

The  thermal  conductivity  of  a  porous  solid  is  not  simply  the  thermal  conductivity 
of  the  solid  material.  The  effective  thermal  conductivity  of  a  saturated  porous  material 
comprised  of  a  packed  bed  of  spheres  is  given  by 

^  k[(2k,+k,)-2(\-s%k,-kj] 

[2*,+t,+(l-e)(t,-*.)] 

where  is  the  effective  thermal  conductivity  of  the  porous  material,  is  the  thermal 
conductivity  of  the  liquid,  k^  is  the  thermal  conductivity  of  the  solid  material,  and  s  is  the 
porosity  of  the  material  (Chi,  1976:50). 
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2.6  Program  Structure 


The  two  subroutines,  TLIQUID  and  TSOLID,  numerically  evaluate  the  energy 

equations  for  the  liquid  coolant  and  porous  solid.  Since  the  two  energy  equations  are  tied 

,  ,  ,  .  ,  .  h,A'{T  -T.)  • 

through  the  internal  convection  term,  — — the  two  temperature  distributions 


Adx 


have  to  be  solved  iteratively.  The  INPUT  subroutine  reads  in  the  initial  parameters  for 
porosity,  solid  thermal  conductivity,  coolant  channel  pressure,  and  solid  surface  area 
exposed  to  coolant  flow  within  the  solid.  The  INPUT  subroutine  also  specifies  an  initial 
temperature  profile  through  the  porous  wall  for  both  the  liquid  and  solid.  The  TLIQUID 
subroutine  uses  a  backward  differencing  technique  and  the  solid  temperature  profile  to 
calculate  a  new  liquid  temperature  profile  through  the  porous  wall.  The  TSOLID 
subroutine  uses  the  updated  liquid  temperature  profile  and  a  forward  differencing 
technique  to  obtain  a  new  solid  temperature  profile  through  the  wall.  At  each  point  a  new 
time  step  is  calculated  based  upon  the  stability  criterion.  After  calculating  the  two 
profiles,  the  new  temperature  values  were  compared  to  values  from  the  previous  time  step. 
Once  the  difference  was  less  than  a  specified  tolerance,  the  solutions  had  converged. 

After  both  temperature  profiles  converge,  the  OUTPUT  subroutine  is  called  to  write  the 
temperature  profiles  to  a  specified  filename. 
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III.  Model  Validation 


This  chapter  provides  a  discussion  of  the  validation  methods  used  for  each  of  the 
subroutines  within  the  numerical  model.  The  chapter  begins  with  a  brief  description  of  the 
validation  of  property  and  heat  transfer  coefficient  subroutines.  The  next  section  describes 
the  results  of  the  various  test  cases  used  to  validate  the  TSOLID  subroutine.  The  chapter 
ends  with  a  description  of  the  TLIQUID  validation  and  results. 

3.1  Validation  of  Property  and  Heat  Transfer  Coefficient  Subroutines 

The  author  verified  the  H2PROP,  HOTGASPROP,  HCOOLSIDE,  HOTGAS, 
HCOOL,  and  FINDTAW  subroutines  by  hand  calculations  and  comparison  with  property 
data.  All  subroutines  generated  correct  results. 

3.2  TSOLID  Validation 

The  author  used  three  incremental  steps  to  validate  the  TSOLID  subroutine.  In 
validating  the  TSOLID  subroutine,  energy  transfer  to  the  liquid  was  ignored  and  the  solid 
was  treated  with  zero  porosity.  The  first  step  was  to  test  the  model  with  two  fixed 
boundary  temperatures.  The  two  temperatures  were  138.89  K  and  3656.7  K.  The  exact 
temperature  profile  in  a  planar  wall  with  no  internal  heat  generation  is  a  straight  line 
governed  by  Equation  (22)  (Incropera  and  DeWitt,  1990:81). 

(22) 

where  7^(x)  is  the  solid  temperature  as  a  function  of  distance,  x,  through  the  wall,  is 
the  cold  wall  temperature,  7^  is  the  hot  wall  temperature,  and  L  is  the  total  thickness  of 
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Temperature  (K) 


the  wall.  The  model  generated  the  profile  shown  in  Figure  2.  As  the  reader  can  see,  the 
solid  temperature  profile  is  linear,  beginning  and  ending  at  the  appropriate  boundary 
temperatures. 


Figure  2.  Fixed  Temperature  Boundary  Condition  Solid  Temperature  Profile 


The  second  step  was  to  impose  two  fixed  convection  boundary  conditions.  The 
adiabatic  wall  temperature,  was  fixed  at  3656.7  K,  the  hot  side  heat  transfer 
coefficient,  was  fixed  at  20,000  W/(m2*K),  the  cold  side  heat  transfer  coefficient, 

was  fixed  at  10,000  W/(m2*K),  and  the  bulk  temperature  of  the  coolant,  7^ ,  was  held 
constant  at  138.89  K.  The  exact  temperature  profile  is  a  straight  line  governed  by 
Equation  (22).  In  this  case,  the  temperatures  7^  and  are  unknown.  Using  the  heat 
flux  correlations  shown  in  Sections  2.2  and  2.3,  and  knowing  that  the  heat  flux  into  the 
wall  equals  the  heat  flux  out  of  the  wall,  the  reader  can  see  that  7^  and  7^^  are  functions 
of  and  7^ .  The  hot  wall  and  cold  wall  temperatures  are  calculated  using  the  equations 
below  (Incropera  and  DeWitt,  1990:83). 


•■WH 


(23) 


(24) 


T 

^wc 


(25) 


The  TSOLED  model  generated  a  solid  temperature  profile  that  matched  the  exact  solution. 
Some  select  temperature  values  are  listed  in  Table  1 . 
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TABLE  1 

FIXED  CONVECTIVE  HEAT  TRANSFER  COEFFICIENT  TEST  CASE  RESULTS 


X  Location 

Exact  Values 

Computed  Values 

millimeters 

K 

K 

0 

2452.752 

2452.377 

0.305816 

2472.966 

2472.599 

0.7112 

2499.774 

2499.412 

The  third  step  was  to  run  the  model  with  realistic  heat  transfer  coefficients 
calculated  from  the  same  material,  combustion,  and  coolant  properties  as  the  SSME  MCC. 
The  model  results  were  compared  to  the  actual  measured  SSME  MCC  temperature 
values.  The  measured  values  for  the  SSME  MCC  wall  temperatures  were  838.7  K  on  the 
hot  gas  side  and  505.4  K  on  the  cold  gas  side  (Cook  and  Coffey,  1973:8).  These  results 
are  shown  in  Figure  3 .  The  dashed  lines  represent  the  maximum  error  caused  by  the  heat 
transfer  correlations  utilizing  the  ±25%  error  on  the  cold  gas  side  and  the  ±2%  error  on 
the  hot  gas  side.  As  the  reader  can  see,  considering  the  possible  error,  the  model 
generates  values  very  close  to  experimental  results.  Based  upon  these  results,  the  author 
is  confident  the  boundary  condition  heat  transfer  equations  and  the  TSOLID  model 
generate  valid  results. 
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Figure  3.  Regeneratively  Cooled  Model  vs.  Experimental  Results 


3.3  TLIOUID  VaUdation 


The  exact  solution  for  constant  liquid  property  flow  through  a  packed  bed  of 


spheres  held  at  a  constant  wall  temperature  is 


5_^  =  exp(-^) 

Ts  -  T,j„  m,c 
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where  7^  is  the  solid  temperature,  T,-„  is  the  inlet  bulk  temperature,  is  the  exit  bulk 

temperature,  h^,  is  the  average  heat  transfer  coefficient,  W;  is  the  mass  flow  rate  of  coolant 
through  the  bed.  A*  is  the  solid  surface  area  exposed  to  the  fluid,  and  is  the  specific 

heat  of  the  fluid  (Incropera  and  DeWitt,  1992:440).  Various  grid  sizes  were  used  to  plot 
the  computed  TLIQUID  temperature  profile  for  a  constant  wall  temperature  of  2000  K 
and  constant  hydrogen  properties  at  138.89  K  and  12.0  MPa.  The  liquid  temperature 
profiles  for  each  of  the  grid  sizes  are  plotted  in  Figure  4. 
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IV.  Results 

This  chapter  summarizes  the  results  for  the  variation  of  parameters  study.  Each 
section  addresses  the  varying  parameter's  effect  on  maximum  wall  temperature  and 
maximum  temperature  gradient  and  compares  these  values  to  the  regeneratively  cooled 
wall  case.  The  maximum  wall  temperature  and  maximum  temperature  gradient  for  the 
regeneratively  cooled  wall  is  838.7  K  and  4.26x10^  K/m  respectively.  To  eliminate 
pressure  differential  dependence,  all  trends  are  non-dimensionalized  by  using  the  blowing 
ratio  of  the  coolant  flow.  The  values  used  in  this  parameter  study  are  presented  below  in 
Table  2. 

TABLE  2 

POROSITIES,  THERMAL  CONDUCTIVITIES,  AND  SPHERE  RADII  USED  IN 

PARAMETER  STLODY. 


Porosity/CO 

0.1546/2.5360 

0.1 

10 

0.3493/1.9519 

1 

50 

0.3686/1.8940 

10 

100 

0.4250/1.7248 

100 

0.4764/1.5708 

1000 
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The  first  section  addresses  some  general  results  pertaining  to  flow  through  a 
porous  medium.  The  second  section  addresses  the  effects  of  varying  thermal  conductivity 
and  blowing  ratio  on  maximum  wall  temperature  and  maximum  temperature  gradient.  The 
third  section  addresses  the  trends  caused  by  varying  solid  porosity  and  blowing  ratio. 
Finally  the  fourth  section  provides  the  effects  varying  porous  sphere  radius  has  upon 
maximum  wall  temperature  and  maximum  temperature  gradient. 

41  General  Result 

A  general  result  can  be  drawn  before  looking  at  the  variation  of  parameters.  In  the 
past,  whenever  scientists  have  considered  flow  through  a  porous  medium,  they  have 
assumed  that  the  liquid  and  solid  are  very  close  to  the  same  temperature  throughout 
(Schneider,  1955:218).  This  study  found  this  to  be  true  for  low  porosity  solids,  but  not 
for  high  porosity  solids.  The  liquid  temperature  remained  within  0. 1%  of  the  solid 
temperature  profile  through  the  entire  wall  thickness  at  a  low  porosity  of  0. 1 546.  This 
case  is  shown  in  Figure  5.  Because  of  the  low  porosity,  the  porous  wall's  effective  thermal 
conductivity  is  highest.  With  a  higher  thermal  conductivity,  the  heat  is  conducted  away 
from  the  wall  surface  and  is  distributed  through  the  wall,  thereby  decreasing  the  maximum 
temperature  gradient. 

At  higher  porosities,  the  effective  thermal  conductivity  is  decreased.  Therefore, 
very  little  heat  is  conducted  through  the  solid  and  the  cooling  effect  of  the  liquid 
dominates.  This  causes  the  solid  and  liquid  to  remain  at  essentially  the  same  temperature 
throughout  most  of  the  wall.  The  two  temperatures  diverge  within  the  last  9%  of  the  wall 
thickness  with  a  maximum  difference  of  39%.  This  effect  is  shown  in  Figure  6. 

The  author  concludes  the  assumption  of  solid  and  liquid  temperatures  nearly  the 
same  throughout  the  porous  wall  is  valid.  This  assumption  breaks  down  in  high  porosity 
solids  near  the  hot  gas  side  boundary. 
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Figure  5.  0.1546  Porous  Solid  and  Liquid  Temperature  Profiles 
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Figure  6.  0.425  Porous  Solid  and  Liquid  Temperature  Profiles 


One  observation  from  the  test  data  reveals  that  varying  thermal  conductivity  had 


the  greatest  effect  in  low  porosity  solids.  Increasing  the  thermal  conductivity  from  0. 1 
W/(m*K)  to  1000  W/(m*K)  in  a  solid  with  a  porosity  of  0. 1546  yielded  a  42%  decrease  in 
maximum  wall  temperature  and  a  95%  decrease  in  maximum  temperature  gradient.  The 


same  increase  in  thermal  conductivity  in  a  solid  with  a  porosity  of  0.4764  yielded  only  a 
3 1%  decrease  in  maximum  wall  temperature  and  a  72%  decrease  in  maximum  temperature 
gradient.  The  effective  thermal  conductivity  is  a  very  strong  function  of  porosity.  As 
porosity  increases,  the  effective  thermal  conductivity  approaches  the  thermal  conductivity 
of  the  liquid.  Therefore,  at  low  porosities  the  solid  thermal  conductivity  will  dominate. 

At  high  porosities  the  liquid  thermal  conductivity  will  dominate.  Therefore,  changes  in 
solid  thermal  conductivity  will  have  a  pronounced  effect  in  low  porosity  solids  and  little 
effect  in  high  porosity  solids.  These  results  are  shown  in  Figures  7  and  8. 
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Figure  7.  0. 1546  Porous  Solid  Temperature  Profile  with  Varying  Thermal  Conductivity 
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Figure  8.  0.4764  Porous  Solid  Temperature  Profile  with  Varying  Thermal  Conductivity 

The  results  of  varying  thermal  conductivity  in  a  porous  solid  with  a  porosity  of 
0.3493,  and  a  sphere  size  of  10  micrometers,  on  maximum  wall  temperature  are  presented 
in  Figure  9. 
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Figure  9.  Varying  Thermal  Conductivity  Effect  on  Maximum  Wall  Temperature 

One  observed  trend  was  a  consistent  decrease  in  maximum  wall  temperature  with 
respect  to  blowing  ratio.  For  a  thermal  conductivity  of  0. 1  W/(m*K),  there  was  a  55% 
decrease  in  maximum  wall  temperature  as  blowing  ratio  went  from  0.06%  to  0.6%.  This 
decrease  is  due  to  two  factors.  One  factor  is  the  increased  advection  of  energy  by  the 
transpiration  coolant.  Another  is  the  decrease  in  heat  flux  caused  by  a  decrease  in  heat 
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transfer  coefficient.  To  illustrate  the  effect  of  blowing  upon  the  maximum  wall 
temperature,  the  author  plotted  two  solid  temperature  profiles  in  Figure  10. 


Figure  10.  Blowing  Effect  on  Liquid  Temperature  Profile 
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One  profile  was  calculated  using  a  heat  transfer  coefficient  neglecting  blowing  effects. 

The  other  profile  was  calculated  using  Chen's  correlation  to  include  blowing  effects.  The 
blowing  effect  alone  caused  a  16%  decrease  in  maximum  wall  temperature. 

Another  noticeable  effect  was  the  decrease  in  effect  of  blowing  ratio  on  maximum 
wall  temperature  with  increasing  thermal  conductivity.  There  was  a  27%  decrease  in 
maximum  wall  temperature  for  the  1000  W/(m*K)  thermal  conductivity  case,  compared  to 
the  55%  decrease  for  the  0.1  W/(m*K)  case.  This  effect  is  related  to  the  increase  in 
effective  thermal  conductivity. 

Additionally,  there  is  only  a  6%  difference  in  maximum  wall  temperatures  for 
thermal  conductivities  between  10  and  1000  W/(m*K),  and  a  1%  difference  between  100 
and  1000  W/(m*K)  solids.  Since  most  modem  materials  fall  within  this  range  of  thermal 
conductivities,  porous  solids  can  be  made  of  existing  materials. 

The  transpiration  cooled  wall  could  achieve  a  lower  maximum  wall  temperature 
than  the  regeneratively  cooled  wall  at  a  blowing  ratio  as  low  as  0.4%  for  thermal 
conductivies  between  10  and  1000  W/(m*K).  At  the  maximum  thermal  conductivity  and 
maximum  blowing  ratio,  the  transpiration  cooled  wall  was  35%  cooler  than  the 
regeneratively  cooled  wall.  Therefore  a  porous  wall  made  of  existing  materials  with  the 
given  properties  could  operate  at  lower  maximum  wall  temperatures  than  the 
regeneratively  cooled  wall. 

The  results  of  varying  thermal  conductivity  in  a  porous  solid  with  a  porosity  of 
0.3493,  and  a  sphere  size  of  10  micrometers,  on  maximum  temperature  gradient  are 
presented  in  Figure  1 1 .  One  observed  trend  was  a  consistent  decrease  in  temperature 
gradient  with  increasing  thermal  conductivity.  There  was  a  73%  decrease  in  maximum 
temperature  gradient  as  thermal  conductivity  went  from  0. 1  to  1000  W/(m*K).  The 
decrease  in  temperature  gradient  is  again  due  to  the  increased  effective  thermal 
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conductivity.  The  increased  effective  thermal  conductivity  allows  the  porous  solid  to 
conduct  the  heat  further  away  from  the  hot  gas  surface.  Unfortunately,  the  lowest 
temperature  gradient  achieved  with  transpiration  cooling  is  still  72  times  greater  than  the 
temperature  gradient  of  regenerative  cooling.  A  possible  tradeoff  would  be  that  the 
decreased  temperatures  would  allow  the  porous  material  to  endure  the  steeper 
temperature  gradient. 
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Figure  1 1 .  Varying  Thermal  Conductivity  Effect  on  Maximum  Temperature  Gradient 
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Another  trend  was  an  increase  in  maximum  temperature  gradient  at  low  blowing 
ratios  and  a  decrease  in  maximum  temperature  gradient  at  high  blowing  ratios.  This  effect 
is  possibly  due  to  two  competing  factors.  As  mentioned  earlier,  as  blowing  ratio 
increases,  the  hot  gas  side  heat  transfer  coefficient  decreases.  This  would  cause  a 
decrease  in  temperature  gradient.  But  an  increase  in  blowing  ratio  causes  the  majority  of 
the  heat  transfer  between  the  solid  and  liquid  to  take  place  closer  to  the  hot  gas  surface 
thereby  increasing  the  temperature  gradient.  Since  existing  materials  fall  in  the  mid  to  high 
thermal  conductivity  range  where  these  temperature  gradient  effects  are  less  pronounced, 
this  trend  is  not  important 

From  these  results,  the  optimum  porous  material  would  be  made  of  the  highest 
thermal  conductivity.  The  material  would  also  operate  at  the  highest  blowing  ratio  to 
obtain  the  lowest  maximum  wall  temperature. 

4.3  Varying  Solid  Porosity  Results 

The  results  of  varying  porosity  in  a  porous  solid  with  a  thermal  conductivity  of  100 
W/(m*K),  and  a  sphere  size  of  10  micrometers,  on  maximum  wall  temperature  are 
presented  in  Figure  12. 
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Figure  12.  Varying  Porosity  Effect  on  Maximum  Wall  Temperature 


A  trend  observed  in  the  previous  section  was  that  increasing  porosity  affected  the  effective 
thermal  conductivity  of  the  porous  solid.  Following  this  reasoning,  it  would  appear  that 
the  porous  material  with  the  lowest  porosity  and  highest  thermal  conductivity  would  have 
the  lowest  maximum  wall  temperature.  But  as  the  reader  can  see  in  Figure  12,  maximum 
wall  temperature  decreases  with  increasing  porosity,  then  increases  again.  The  minimum 
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wall  temperature  occurs  near  a  porosity  of  0.3493,  Apparently  there  is  an  inflection  point 
where  the  increased  coolant  available  for  energy  absorption  at  larger  porosities  is  offset  by 
decreased  effective  thermal  conductivity  and  decreased  porous  solid  surface  area.  As 
porosity  increases,  the  surface  area  exposed  to  flow  decreases.  The  optimum  case  seems 
to  be  where  the  porosity  is  large  enough  to  allow  a  large  mass  transfer  of  coolant,  but  not 
large  enough  to  diminish  the  effective  thermal  conductivity.  The  porous  solid  with  a 
porosity  of  0.3493  operates  at  a  15%  lower  temperature  than  the  0.1546  solid  and  at  a 
23%  lower  temperature  than  the  0.4764  solid. 

The  results  of  varying  porosity  in  a  porous  solid  with  a  thermal  conductivity  of  100 
W/(m*K),  and  a  sphere  size  of  10  micrometers,  on  maximum  temperature  gradient  are 
presented  in  Figure  13.  As  porosity  increases,  the  maximum  temperature  gradient 
increases.  Since  the  effective  thermal  conductivity  decreases  with  increasing  porosity,  the 
amount  of  heat  conducted  into  the  porous  solid  decreases.  Therefore  all  the  heat  transfer 
between  the  transpiration  coolant  and  the  solid  has  to  occur  near  the  hot  gas  side  wall. 
Therefore  the  temperature  gradients  have  to  increase  because  the  temperature  change  has 
to  occur  over  a  shorter  and  shorter  distance.  The  temperature  gradient  increased  over  27 
times  as  porosity  increased  from  0. 1546  to  0.4764. 

From  these  results,  if  only  maximum  wall  temperature  is  considered,  the  optimum 
porous  material  would  be  made  with  a  porosity  of  approximately  0.3493.  The  design 
engineer  would  have  to  consider  the  effect  of  increased  temperature  gradients  due  to  the 
higher  porosity. 
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Figure  13.  Varying  Porosity  Effect  on  Maximum  Temperature  Gradient 


4  4  Varying  Porous  Sphere  Radius  Results 

The  porous  material  sphere  sizes  used  were  limited  due  to  the  wall  thickness. 
Because  the  wall  is  only  71 1.2  micrometers  thick,  spheres  larger  than  200  micrometers  in 
diameter  were  not  considered.  The  results  of  varying  sphere  size  in  a  porous  solid  with  a 


40 


Maximum  Wall  Temperature  (K) 


porosity  of  0.3493,  and  a  thermal  conductivity  of  100  W/(m*K),  on  maximum  wall 
temperature  are  presented  in  Figure  14. 
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Figure  14.  Varying  Porous  Sphere  Radius  Effect  on  Maximum  Wall  Temperature 


There  was  a  142%  increase  in  maximum  wall  temperature  by  increasing  the  sphere  radius 
from  10  micrometers  to  100  micrometers.  This  effect  is  due  to  the  drastically  reduced 
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surface  area  exposed  to  the  transpiration  coolant  flow.  Since  the  surface  area  exposed  to 
the  coolant  flow  takes  the  form 


A'  =  CO 


\^pJ 


(27) 


where  A"  is  the  solid  surface  area  exposed  to  the  flow,  CO  is  a  proportionality  constant 
defined  by  the  geometrical  arrangement  of  the  spherical  particles,  x  is  the  distance 
through  the  wall,  is  the  spherical  particle  radius,  and  A  is  the  transpiration  coolant  flow 
cross-sectional  area.  From  this  equation.  A*  increases  as  decreases.  Therefore,  the 

most  heat  transfer  fi-om  the  solid  to  the  transpiration  coolant  occurs  with  the  smallest 
radius  spheres. 

The  results  of  varying  sphere  size  in  a  porous  solid  with  a  porosity  of  0.3493,  and 
a  thermal  conductivity  of  100  W/(m*K),  on  maximum  temperature  gradient  are  presented 
in  Figure  15. 


42 


Maximum  Temperature  Gradient  (K/m) 


4.00E+7  — 

O.OOE+0  — - 

0,00E+0 


2,00E-3 


4.00E-3  6.00E-3 

Blowing  Ratio 


8.00E-3  1.00E-2 


Figure  15.  Varying  Porous  Sphere  Radius  Effect  on  Maximum  Temperature  Gradient 


Maximum  temperature  gradient  varied  less  than  10%  across  the  entire  range  of  blowing 
ratios.  Therefore,  there  was  little  to  no  effect  of  varying  porous  sphere  radius  on 
temperature  gradient.  From  these  results,  the  optimum  porous  material  would  be  made  of 
the  smallest  diameter  spheres  possible  Avithout  impeding  coolant  flow. 
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V.  Conclusions  and  Recommendations 


5.1  Conclusions 

The  author  created  a  numerical  model  to  compare  the  effect  of  porosity,  thermal 
conductivity,  and  blowing  ratio  upon  the  solid  temperature  profile  within  a  porous, 
transpiration  cooled  rocket  nozzle.  The  maximum  wall  temperature  and  maximum 
temperature  gradient  for  each  case  was  compared  to  the  maximum  wall  temperature  and 
maximum  temperature  gradient  for  the  regeneratively  cooled  wall  case.  The  transpiration 
cooled  wall  could  achieve  lower  maximum  wall  temperatures  than  the  regeneratively 
cooled  wall  with  blowing  ratios  as  low  as  0.4%.  The  transpiration  cooled  wall  could 
operate  at  a  maximum  35%  colder  than  the  regeneratively  cooled  wall.  Unfortunately, 
transpiration  cooling  results  in  a  temperature  gradient  that  is  72  times  steeper  than  the 
regeneratively  cooled  wall.  Therefore,  the  material  chosen  for  the  porous  wall  must  be 
able  to  withstand  extreme  temperature  gradients.  The  maximum  temperature  gradient  is 
located  at  the  hot  gas  surface  of  the  porous  wall.  The  optimum  configuration  for  the  cases 
studied  would  be  a  medium  (0.3493)  porosity,  high  (1000  W/(m*K))  thermal 
conductivity,  and  high  (1%)  blowing  ratio  transpiration  cooled  nozzle.  The  thermal 
conductivity  and  blowing  ratio  values  represent  the  largest  values  tested  in  this  study. 

Since  the  medium  thermal  conductivity  results  were  not  substantially  different  from  the 
high  thermal  conductivity  results,  existing  materials  (100  -  500  W/m*K  thermal 
conductivities)  could  be  used  to  make  the  porous  nozzle  wall. 
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5  2  Recommendations 

The  author  hopes  that  this  study  will  shed  some  new  light  on  the  feasibility  of 
using  transpiration  cooling  as  a  method  for  cooling  the  throat  region  of  long  operating 
duration,  liquid-fueled  rocket  nozzles.  One  recommendation  is  to  build  upon  the  model 
shown  in  this  study  and  tailor  it  to  other  specific  problems.  Modifications  could  include 
increasing  the  accuracy  of  the  liquid  initial  boundary  condition  by  solving  the  turbulent 
boundary  layer  problem  to  obtain  the  two  coolant  temperature  values  prior  to  contact  with 
the  porous  solid.  Another  recommedation  would  be  to  choose  more  accurate  heat  transfer 
correlations  for  the  boundary  conditions.  This  thesis  test  case  using  the  SSME  MCC  is 
specifically  narrow.  The  blowing  ratios  used  were  small,  and  geometric  and 
thermodynamic  conditions  were  tailored  for  the  SSME  MCC  case.  The  author 
recommends  that  future  research  could  use  the  basic  outline  of  this  model  to  test  against 
experimental  transpiration  cooled  thrust  chambers. 
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Appendix  A.  Derivation  of  Energy  Equations 


This  appendix  presents  the  derivation  of  the  energy  equations  used  to  calculate  the 
solid  and  liquid  temperature  profiles  and  their  conversion  to  the  finite  difference  forms. 

The  first  section  addresses  the  solid  energy  equation.  The  second  section  covers  the  liquid 
energy  equation.  The  third  and  fourth  sections  address  the  finite  difference  forms  of  the 
solid  and  liquid  energy  equations  respectively. 

A.  1  Solid  Energy  Equation 

The  derivation  begins  by  performing  an  energy  balance  on  a  differential  control 
volume  of  the  porous  solid  material.  The  energy  balance  takes  the  form 

^in  +  ^generated  ~  ^out  ~  ^stored  (Al) 

where  is  the  rate  of  energy  entering  the  control  volume,  the  rate  of  energy 

generated  within  the  control  volume,  is  the  rate  of  energy  leaving  the  control  volume, 
and  ‘s  t^e  rate  of  energy  stored  within  the  control  volume  (Incropera  and  DeWitt, 
1990;  15).  This  study  only  considers  heat  transfer  in  one  dimension  represented  by  x,  with 
no  heat  generation  within  the  solid.  The  energy  processes  within  the  porous  solid  are 
represented  below  in  Figure  Al. 
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x+dx 


Figure  Al.  Energy  Balance  on  a  Solid  Differential  Control  Volume 


Representing  the  energy  processes  in  equation  form  below  yields  the  equation  below. 


-k,  —  A-  -k^^A  +  ~{-k^^^Acbc  -h.A\T-T,)^—(pcT)Adx  (A2) 

where  k^  is  the  effective  thermal  conductivity  of  the  porous  solid  wall,  — ^  is  the 

dx 

temperature  gradient  through  the  porous  solid  wall,  A  is  the  conduction  cross  sectional 
area,  cbc  is  the  differential  distance  in  the  x  direction,  h^,  is  the  heat  transfer  coefficient 
between  the  porous  solid  and  the  transpiration  coolant.  A'  is  the  surface  area  of  the 
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transient  temperature  effects  within  the  porous  solid. 

Since  the  energy  transfer  is  one  dimensional  in  space  and  time,  all  partial 
differentials  become  total  derivatives.  Also,  since  the  effective  thermal  conductivity  takes 
into  account  the  porosity  of  the  solid  material,  the  conduction  cross  sectional  area  used 
becomes  the  total  area.  A,  not  the  (l  -  £)y4  area  as  previous  analyses  have  assumed 

(Schneider,  1955:219).  Assuming  constant  properties  with  respect  to  the  x  direction  and 
time,  canceling  like  terms,  and  dividing  by  Adx  yields  the  final  form  of  the  solid  energy 
equation  shown  below. 


Adx  ^  ’’’  dt 


(A3) 


A.  2  Liquid  Energy  Equation 

The  derivation  of  the  liquid  energy  equation  also  begins  with  an  energy  balance  on 
a  differential  control  volume  in  the  x  direction.  The  liquid  energy  equation  assumes  there 
is  no  heat  generation  within  the  liquid  and  that  the  liquid  behaves  like  an  ideal  gas.  The 
energy  processes  that  take  place  within  the  liquid  inside  the  porous  solid  are  represented  in 
Figure  A2. 
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This  energy  balance  assumes  that  the  only  work  done  on  the  fluid  is  by  pressure  forces  and 
body  forces  are  negligible.  The  area  used  in  these  equations  is  the  liquid  cross  sectional 
area  defined  by  the  porosity /f;  =  sA .  The  advection  term  represents  the  internal  and 
kinetic  energy  of  the  fluid  flow  crossing  the  control  volume.  A  more  convenient  form  uses 
enthalpy  rather  than  internal  energy.  Enthalpy  is  defined  as 


(A4) 
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where  i  is  the  enthalpy  of  the  fluid,  e  is  the  internal  energy  of  the  fluid,  P  is  the  pressure 
on  the  fluid,  and  is  the  density  of  the  fluid  (Incropera  and  DeWitt,  1990:334). 


Combining  terms  from  all  the  energy  processes,  assuming  steady  state,  assuming  the 


u, 


coolant  is  an  ideal  gas,  assuming  the  p,  term  is  small  compared  to  the  temperature 


gradient  within  the  fluid,  and  further  simplifying  yields  the  final  form  of  the  energy 
equation  shown  below. 


‘  dx^  A,dx 


^  =  0 
dx 


(A5) 


A.3  Finite  Difference  Form  of  the  Solid  Energy  Equation 

This  study  employed  an  explicit  finite  differencing  technique  to  solve  the  solid 
differential  energy  equation  for  the  solid  temperature  profile.  Dividing  the  porous  wall 
into  N  regions  means  that  there  has  to  be  N+1  grid  points.  The  first  grid  point  represents 
the  conditions  on  the  cold  gas  side  of  the  wall.  Consequently,  the  N+1  grid  point 
represents  the  hot  gas  side  conditions.  At  each  grid  point,  the  differential  equation 
requires  initial  conditions  in  space  and  time.  Because  the  problem  is  second  order  in  space 
and  first  order  in  time,  the  problem  requires  two  initial  conditions  in  space  and  one  initial 
condition  in  time.  The  general  form  of  the  finite  difference  form  is  given  below  in 
Equation  A6. 


[7r(JV  +  l)+7r“(A'-l)-27r"(W)]  hTA-{Tt‘(N)-V“(N)\ 


dx^ 

old  i 


Adx 


Pf 


s  p, 


[Tr(N)-Tr{N)] 

dt 


(A6) 
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where  the  (A'^  + 1),  (A^  - 1),  and  (N)  terms  refer  to  the  grid  point  being  used  for  the 
calculation,  and  the  old  and  new  superscripts  refer  to  the  time  dependence  of  the  equation 
(Anderson  and  others,  1984:333).  The  reader  can  see  that  the  equation  can  be  solved  for 
directly  as  long  as  the  old  temperatures  at  the  ( A  + 1),  (N  - 1),  and  (N)  are 
known.  This  equation  has  to  be  modified  at  the  first  and  last  grid  points  because  there  is 
not  two  sided  conduction  taking  place.  The  boundary  conditions  at  these  locations  are 
governed  by  the  heat  flux  between  the  porous  solid  and  the  hot  combustion  gases  and  the 
cold  convective  coolant  respectively.  The  TSOLED  subroutine  begins  with  an  initial 
temperature  distribution  at  time  equals  zero.  The  subroutine  selects  a  time  step,  dt,  based 
upon  the  grid  size  and  stability  criteria,  and  marches  forward  through  the  porous  wall  in 
both  space  and  time  until  the  temperature  does  not  change  beyond  a  threshold  value. 

Once  this  tolerance  is  reached,  the  solid  has  reached  its  steady  state  temperature  profile. 
A.4  Finite  Difference  Form  of  the  Liquid  Energy  Equation 

Unfortunately,  the  liquid  solution  is  a  bit  more  complicated.  Unlike  the  solid 
energy  equation,  the  second  spatial  initial  condition  is  not  known  at  the  last  grid  point. 
Fortunately,  two  coolant  temperatures  are  known  just  before  entering  the  porous  solid. 
This  study  assumes  these  temperatures  to  be  the  bulk  fluid  temperature  of  the  coolant 
within  the  coolant  channel.  More  accurate  values  could  be  found  by  calculating  the 
boundary  layer  temperature  profile  within  the  coolant  channel.  For  the  purposes  of  this 
study,  the  bulk  temperature  of  the  fluid  was  sufficient.  The  finite  difference  form  of  the 
liquid  energy  equation  takes  a  slightly  modified  form  from  the  solid  energy  equation. 


,  {Tr{M-2)  +  Tr{M)-2Tr'{M-r)\  , 

ICj  -  I 

dx  A;dx 


(A7) 


Pi 


dx 


=  0 
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where  the  liquid  values  are  calculated  along  M  points.  N  and  M  are  related  by  M=N+2. 
Therefore,  at  the  first  solid  grid  point  N=l,  the  liquid  temperature  profile  is  M=3.  This 
equation  can  be  solved  directly  for  since  the  values  of  -  2)  and 

1)  are  known  at  all  times.  The  TLIQUID  subroutine  takes  an  initial  solid 
temperature  profile  and  the  two  known  grid  points,  and  and  marches 

forward  through  the  porous  solid  until  it  reaches  the  last  grid  point.  Like  the  TSOLID 
subroutine,  the  TLIQUID  subroutine  compares  the  new  liquid  temperature  profile  to  the 
previous  profile  and  watches  for  the  difference  to  be  within  a  specified  tolerance.  The 
liquid  has  reached  its  steady  state  temperature  profile  once  this  tolerance  is  reached. 
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Appendix  B.  Thesis  Computer  Code 


Attached  is  the  final  computer  code  in  FORTRAN  used  to  compute  the 
temperature  profiles  for  the  porous  solid  and  liquid  transpiration  coolant. 


53 


PROGRAM  THESIS 


C 

C  THIS  PROGRAM  CALCULATES  THE  TEMPERATURE  PROFILES 
C  FOR  A  POROUS  SOLID  AND  THE  LIQUID  TRANSPIRATION 
C  COOLANT  OF  A  ROCKET  NOZZLE  WITH  THE  SPACE  SHUTTLE 
C  MAIN  ENGINE  MAIN  COMBUSTION  CHAMBER  COMBUSTION  PRODUCTS, 

C  HYDROGEN  COOLANT.  AND  NOZZLE  GEOMETRY 
C 

IMPLICIT  NONE 
C 

C  VARIABLES  USED  IN  MAIN  PROGRAM 
C 

C  TSOLD  =  VECTOR  OF  POROUS  SOLID  TEMPERATURE  AT  PREVIOUS  TIME  STEP 
C  TLOLD  =  VECTOR  OF  TRANSPIRATION  COOLANT  TEMPERATURE  AT  PREVIOUS 
C  TIME  STEP 

C  TSNEW  =  SOLID  TEMPERATURE  AT  NEW  TIME  STEP 
C  TLNEW  =  COOLANT  TEMPERATURE  AT  NEW  TIME  STEP 
C  DX  =  X  DISTANCE  STEP  SIZE  THROUGH  WALL 
C  RHOLAVG  =  AVERAGE  COOLANT  DENSITY 
C  MULAVG  =  AVERAGE  COOLANT  VISCOSITY 
C  DIFF  =  TOLERANCE  USED  TO  DETERMINE  CONVERGENCE 
C  DELTAS  =  SUM  OF  THE  SQUARES  DIFFERENCE  OF  SOLID  TEMPERATURES 
C  DELTAL  =  SUM  OF  THE  SQUARES  DIFFERENCE  OF  COOLANT  TEMPERATURES 
C  KL  =  LOCAL  COOLANT  THERMAL  CONDUCTIVITY 
C  DIFFLIM  =  SPECIFIED  TOLERANCE  USED  TO  DETERMINE  CONVERGENCE 
C  N  =  INDEX  USED  FOR  SOLID  TEMPERATURE  ARRAY 
C  M  =  INDEX  USED  FOR  COOLANT  TEMPERATURE  ARRAY 
C  LASTX  =  GRID  SIZE 

C  STEP  =  COUNTER  THAT  INDICATES  NUMBER  OF  ITERATIONS  OF  MAIN  LOOP 
C 

C  DECLARE  THE  VARIABLES 
C 

DOUBLE  PRECISION  TSOLD,  TLOLD,  TSNEW,  TLNEW,  DX 
DOUBLE  PRECISION  RHOLAVG,  MULAVG,  DIFF 
DOUBLE  PRECISION  DELTAS,  DELTAL,  KL,  DIFFLIM 

INTEGER  N,  M,  LASTX,  STEP 
C 

C  DIMENSION  THE  VARIABLES  THAT  ARE  ARRAYS 
C 

DIMENSION  TSOLD(20000),  TLOLD(20000),  TSNEW(20000),  TLNEW(20000) 
DIMENSION  KL(20000) 

C 

C  THIS  COMMON  BLOCK  CONTAINS  THE  GRID  DIMENSIONS  USED 
C 

COMMON  /DIMEN/  LASTX,  DX 
C 

C  BEGIN  THE  PROGRAM  BY  SPECIFYING  THE  INITIAL  CONDITIONS 
C 

CALL  INPUT(TSOLD,  TLOLD,  RHOLAVG,  MULAVG,  DIFFLIM,  TLNEW) 
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INITIALIZE  THE  CONVERGENCE  CRITERIA  VARIABLE,  DIFF.  TO  A  VERY 
LARGE  NUMBER  TO  BEGIN  MAIN  LOOP 

DIFF  =  9999.0 
STEP  =  0 


BEGIN  THE  MAIN  LOOP  AND  CONTINUE  TO  ITERATE  UNTIL  THE  SUM  OF 
THE  SQUARES  OF  THE  SOLID  AND  COOLANT  TEMPERATURES  IS  LESS  THAN 
THE  SPECIFIED  TOLERANCE,  DIFFLIM 

DO  WHILE  (DIFF.  GE.  (DIFFLIM)) 

CALCULATE  A  NEW  COOLANT  TEMPERATURE  PROFILE  BASED  ON  PREVIOUS 
TIME  STEP  SOLID  TEMPERATURE  PROFILE.  ALSO  GENERATE  NEW  COOLANT 
THERMAL  CONDUCTIVITIES,  AVERAGE  VISCOSITY,  AND  AVERAGE  DENSITY 

CALL  TLIQUID(TSOLD,  TLOLD,  KL,  MULAVG,  RHOLAVG,  TLNEW) 

CALCULATE  A  NEW  SOLID  TEMPERATURE  PROFILE  BASED  ON  NEW  COOLANT 
TEMPERATURE  PROFILE,  AND  COOLANT  THERMAL  CONDUCTIVITIES 

CALL  TSOLID(TSOLD,  TLOLD,  KL,  TSNEW) 

INITIALIZE  THE  SUM  OF  THE  SQUARES  DIFFERENCES,  DELTAS  AND  DELTAL 
TO  ZERO  BEFORE  SUMMING  TERMS 

DELTAS  =  0.0 
DELTAL  =  0.0 

SUM  UP  THE  SQUARES  OF  THE  DIFFERENCES  BETWEEN  THE  SOLID  AND  COOLANT 
TEMPERATURES 

DON=I,(LASTX+I) 

M  =  N  +  2 

DELTAS  =  DELTAS  +  ((TSNEW(N)-TSOLD(N))**2.0) 

DELTAL  =  DELTAL  +  ((TLNEW(M)-TLOLD(M))**2.0) 

END  DO 

THE  TOLERANCE  VARIABLE,  DIFF,  IS  ASSIGNED  THE  SUM  OF  THE  SQUARES 
DIFF  =  DELTAS  +  DELTAL 

ASSIGN  THE  NEW  TIME  STEP  VALUES  TO  THE  PREVIOUS  TIME  STEP  VARIABLES 

DON=L(LASTX+I) 

M  =  N  +  2 

TSOLD(N)  =  TSNEW(N) 

TLOLD(M)  =  TLNEW(M) 

END  DO 
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C  INCREMENT  THE  ITERATION  COUNTER 
C 

STEP  =  STEP+I 
C 

C  END  THE  MAIN  LOOP 
C 

END  DO 
C 

C  SEND  THE  NEWEST  TEMPERATURE  PROFILES  TO  THE  OUTPUT  SUBROUTINE  TO 
C  BE  STORED  IN  A  FILE 
C 

CALL  OUTPUT(TSNEW,  TLNEW) 

C 

C  END  THE  MAIN  PROGRAM 
C 

STOP 

END 

C 

C  THIS  SUBROUTINE  PROVIDES  THE  INITIAL  CONDITIONS  AND  CONSTANTS 
C  USED  IN  THE  REST  OF  THE  PROGRAM 
C 

SUBROUTINE  INPUT(TSOLD,  TLOLD,  RHOLAVG,  MULAVG,  DIFFLIM,  TLNEW) 

IMPLICIT  NONE 
C 

C  VARIABLES  USED  IN  INPUT 
C 

C  TSOLD  =  VECTOR  OF  POROUS  SOLED  TEMPERATURE  AT  PREVIOUS  TIME  STEP 
C  TLOLD  =  VECTOR  OF  TRANSPIRATION  COOLANT  TEMPERATURE  AT  PREVIOUS 
C  TIME  STEP 

C  EPS  =  POROSITY  OF  THE  POROUS  SOLID 
C  R  =  POROUS  SPHERE  RADIUS 

C  TB  =  BULK  TEMPERATURE  OF  THE  COOLANT  IN  THE  COOLANT  CHANNEL 
C  DX  =  X  DISTANCE  STEP  SIZE  THROUGH  WALL 

C  CO  =  PROPORTIONALITY  CONSTANT  USED  TO  DETERMINE  SOLID  SURFACE  AREA 
C  MDOTC  =  MASS  FLOW  RATE  OF  COOLANT  IN  COOLANT  CHANNEL 
C  DH  =  HYDRAULIC  DIAMETER  OF  COOLANT  CHANNEL 
C  PB  =  COOLANT  LINE  PRESSURE  (BULK  CONDITIONS) 

C  PI  =  3.14159 

C  THROATD  =  NOZZLE  AERODYNAMIC  THROAT  DIAMETER 
C  MDOTF  =  MASS  FLOW  RATE  OF  COMBUSTION  PRODUCTS 
C  W  =  VISCOSITY,  TEMPERATURE  EXPONENT 
C  PH  =  HYDRAULIC  PERIMETER  OF  COOLANT  CHANNEL 
C  GAM  =  RATIO  OF  SPECIFIC  HEATS  OF  COMBUSTION  PRODUCTS 
C  TOG  =  STAGNATION  TEMPERATURE  OF  COMBUSTION  PRODUCTS 
C  TG  =  STATIC  TEMPERATURE  OF  COMBUSTION  PRODUCTS  AT  THROAT 
C  DELTAP  =  CHANGE  IN  PRESSURE  ACROSS  THE  POROUS  WALL 
C  THICKNESS  =  POROUS  NOZZLE  WALL  THICKNESS  AT  THROAT 
C  PERMEABLE  =  POROUS  WALL  PERMEABILITY 
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C  RHOLAVG  =  AVERAGE  COOLANT  DENSITY 
C  MULAVG  =  AVERAGE  COOLANT  VISCOSITY 
C  PLAVG  =  AVERAGE  PRESSURE  WITHIN  POROUS  WALL 
C  PTHROAT  =  STATIC  THROAT  PRESSURE 
C  RHOB  =  COOLANT  DENSITY  AT  BULK  CONDITIONS 
C  DP  =  PRESSURE  STEP  SIZE  THROUGH  WALL 
C  RHOS  =  SOLID  DENSITY 
C  CPS  =  SOLID  SPECIFIC  HEAT 
C  KW  =  SOLID  THERMAL  CONDUCTIVITY 
C  MUB  =  COOLANT  VISCOSITY  AT  BULK  CONDITIONS 
C  KB  =  COOLANT  THERMAL  CONDUCTIVITY  AT  BULK  CONDITIONS 
C  CPB  =  COOLANT  SPECIFIC  HEAT,  CONSTANT  PRESSURE,  AT  BULK  CONDITIONS 
C  DIFFLIM  =  SPECIFIED  TOLERANCE  USED  TO  DETERMINE  CONVERGENCE 
C  PBB  =  COOLANT  CHANNEL  PRESSURE  FOR  REGENERATIVE  COOLING  CASE 
C  TLAVG  =  AVERAGE  COOLANT  TEMPERATURE  THROUGH  POROUS  WALL 
C  MUOF  =  COMBUSTION  PRODUCTS  VISCOSITY  AT  TOG 
C  CPOF  =  COMBUSTION  PRODUCTS  SPECIFIC  HEAT  AT  TOG 
C  PROF  =  COMBUSTION  PRODUCTS  PRANDTL  NUMBER  AT  TOG 
C  TLNEW  =  COOLANT  TEMPERATURE  AT  NEW  TIME  STEP 
C  KLAVG  =  COOLANT  THERMAL  CONDUCTIVITY  AT  TLAVG  AND  PLAVG 
C  CPLAVG  =  COOLANT  SPECIFIC  HEAT  AT  TLAVG  AND  PLAVG 
C  TLSTART  =  INITIAL  COOLANT  TEMPERATURE  AT  FIRST  GRID  POINT 
C  TLEND  =  INITIAL  COOLANT  TEMPERATURE  AT  LAST  GRID  POINT 
C  DTEMP  =  INCREMENTAL  CHANGE  IN  TEMPERATURE  FOR  INITIAL  PROFILE 
C  HCL  =  PACKED  BED  OF  SPHERES  HEAT  TRANSFER  COEFFICIENT 
C  MDOTL  =  MASS  FLOW  RATE  OF  COOLANT  THROUGH  POROUS  WALL 
C  N  =  INDEX  USED  FOR  SOLED  TEMPERATURE  ARRAY 
C  M  =  INDEX  USED  FOR  COOLANT  TEMPERATURE  ARRAY 
C  LASTX  =  GRID  SIZE 

C  FILE  1  =  FILENAME  USED  TO  STORE  OUTPUT 
C 

C  DECLARE  THE  VARIABLES 
C 

DOUBLE  PRECISION  TSOLD,  TLOLD,  EPS,  R,  TB,  DX,  CO 
DOUBLE  PRECISION  MDOTC,  DH,  PB,  PI,  THROATD,  MDOTF,  W,  PH 
DOUBLE  PRECISION  GAM,  TOG,  TG,  DELTAP,  THICKNESS,  PERMEABLE 
DOUBLE  PRECISION  RHOLAVG,  MULAVG,  PLAVG,  PTHROAT,  RHOB 
DOUBLE  PRECISION  DP,  RHOS,  CPS,  KW,  MUB,  KB,  CPB,  DIFFLIM 
DOUBLE  PRECISION  PBB,  TLAVG,  MUOF,  CPOF,  PROF,  TLNEW 
DOUBLE  PRECISION  KLAVG,  CPLAVG,  TLSTART,  TLEND,  DTEMP 
DOUBLE  PRECISION  HCL,  MDOTL 

INTEGER  LASTX,  N,M 

CHARACTER*64  FILEl 

DIMENSION  THE  VARIABLES  THAT  ARE  ARRAYS 

DIMENSION  TSOLD(20000),  TLOLD(20000),  TLNEW(20000) 

C 
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TfflS  COMMON  BLOCK  CONTAINS  THE  GRID  DIMENSIONS  USED 
COMMON  /DIMEN/  LASTX,  DX 

THIS  COMMON  BLOCK  CONTAINS  THE  OUTPUT  FILENAME 
COMMON  /FILEBLK/  FILE  I 

THIS  COMMON  BLOCK  CONTAINS  VARIABLES  COMMON  TO  TSOLID  AND  TLIQUID 

COMMON  /TSTL/  CO,  EPS,  R,  PB,  TB,  KW,  HCL,  MDOTL 

THIS  COMMON  BLOCK  CONTAINS  VARIABLES  USED  IN  TLIQUID  ONLY 

COMMON  /TLBLK/  PLAVG,  TLAVG,  RHOB,  DP,DELTAP,PERMEABLE,THICKNESS 

THIS  COMMON  BLOCK  CONTAINS  VARIABLES  USED  IN  TSOLID  ONLY 

COMMON  /TSBLK/  RHOS,  CPS,  MDOTC,  DH,  PI,  THROATD,  MDOTF 
COMMON  /TSBLK/  PH,  W,  GAM,  TOG,  TG,  PBB,  MUOF,  CPOF,  PROF 

SPECIFY  CONSTANT  VALUES 

PTHROAT  =  I L890000000D6 

PBB  =  4L37D6 

THICKNESS  =  7. 1 12D-4 

TB=  138.89 

TOG  =  3656.7 

TG  =  3439.4 

W  =  0.875 

GAM  =  1.1481 

MDOTF  =  468.3 134 

THROATD  =  0.26006 

PI  =  3.14159 

DH=  L45143D-03 

PH  =  7.112D-03 

MDOTC  =  3.344D-02 

RHOS  =  19300.0 

CPS  =  129.0 

MU0F=  1.10236699599D-4 
CPOF  =  3740.68492 
PROF  =  0.6151 

READ  IN  VARIABLES  SPECIFIC  TO  THIS  RUN 

READ  (*.1)FILE1 
READ  (*,*)  DIFFLIM 
READ  (*,*)  LASTX 
READ  (*,*)  PB 
READ  (*,*)  R 
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READ  (*,*)  CO 
READ  (*,*)  EPS 
READ  (*,*)  KW 

GENERATE  HYDROGEN  PROPERTIES  AT  BULK  CONDITIONS 
CALL  H2PROP(PB.  TB,  MUB,  KB,  RHOB,  CPB) 

CALCULATE  THE  REST  OF  THE  NECESSARY  VARIABLES  FROM  KNOWN  VALUES 

PLAVG  =  (PB+PTHROAT)/2.0 
DX  =  THICKNESS/LASTX 
DELTAP  =  (PB-PTHROAT) 

DP  =  DELTAP/LASTX 

PERMEABLE  =  ((2.0*R)**2.0)*(EPS**3.0)/(I50.0*((1.0-EPS)**2.0)) 

TLSTART  =  138.89 
TLEND  =  700.0 

DTEMP  =  (TLEND-TLSTART)/LASTX 

SPECIFY  THE  BEGINNING  LIQUID  AND  SOLID  TEMPERATURE  PROFILES 

DON=l,(LASTX+l) 

M  =  N  +  2 

TLOLD(M)  =  TLSTART  +  DTEMP*CN-1) 

TSOLD(N)  =  TLOLD(M)+15 
TLNEW(M)  =  TLOLD(M) 

END  DO 

SPECIFY  NEW  AND  OLD  COOLANT  TEMPERATURES  JUST  PRIOR  TO  THE 
POROUS  WALL  SURFACE 

TLOLD(l)  =  TB 
TLOLD(2)  =  TB 
TLNEW(l)  =  TB 
TLNEW(2)  =  TB 

CALCULATE  THE  AVERAGE  COOLANT  TEMPERATURE  WITHIN  THE  POROUS  WALL 
TLAVG  =  (TLOLD(3)+TLOLD(LASTX+3))/2.0 

GENERATE  HYDROGEN  PROPERTIES  AT  THE  AVERAGE  PRESSURE  AND  TEMPERATURE 

CALL  H2PROP(PLAVG,  TLAVG,  MULAVG,  KLAVG,  RHOLAVG,  CPLAVG) 

FORMAT  STATEMENT  FOR  THE  READ 

FORMAT(2(A)) 

RETURN 

END  INPUT  SUBROUTINE 
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C 

END 

TfflS  SUBROUTINE  GENERATES  HYDROGEN  THERMODYNAMIC  PROPERTIES  AT  A 
GIVEN  PRESSURE  AND  TEMPERATURE 

SUBROUTINE  H2PROP(P,  TL,  MU,  K,  RHO,  CP) 

IMPLICIT  NONE 

VARIABLES  USED  IN  H2PROP 

P  =  HYDROGEN  PRESSURE 
TL  =  LOCAL  HYDROGEN  TEMPERATURE 
MU  =  HYDROGEN  VISCOSITY 
K  =  HYDROGEN  THERMAL  CONDUCTIVITY 
RHO  =  HYDROGEN  DENSITY 

CP  =  HYDROGEN  SPECIFIC  HEAT  AT  CONSTANT  PRESSURE 

DECLARE  VARIABLES 

DOUBLE  PRECISION  P,  TL,  MU,  K,  RHO,  CP 

CURVE  FITS  TO  CALCULATE  PROPERTIES 

MU  =  (I.287I7D-I5)*(TL**3.0)+ 

+  (-8.36499D-12)*(TL**2,0)+(2.6I452D-8*TL)+2.38468D-6 
K  =  (-6.67524D-14)*(TL**4.0)+(2.77669D-10)*(TL**3.0)+ 

+  (-3.227I4D-7)*(TL**2.0)+(0.00056901*TL)+0.0729104 
RHO  =  P/(4157.5*TL) 

CP  =  (-I.39859D-11)*(TL**4.0)+(L48717D-7)*(TL**3.0)- 
+  (0.000849414)*(TL**2.0)+(3.74545*TL)+121 16.7 

THESE  VALUES  REPRESENT  PROPERTIES  AT  600K  AND  1 1.89D6  PA 


MU=  L5338D-5 
K  =  .349464738 
RHO  =  4.76648627 
CP  =  14088.49147502 

THESE  VALUES  REPRESENT  PROPERTIES  AT  138.89K  AND  I2.06D6  PA 


MU=  5.858D-6 
K=  .146434006 
RHO  =  20.885461299 
CP=  12620.913431226 
RETURN 


END  H2PROP  SUBROUTINE 
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C 

C  TfflS  SUBROUTINE  CALCULATES  THE  SOLID  TEMPERATURE  PROFILE  GIVEN 
C  THE  PREVIOUS  TIME  STEP  SOLID  AND  COOLANT  TEMPERATURE  PROFILES.  AND 
C  THE  LIQUID  THERMAL  CONDUCTIVITIES 
C 

SUBROUTINE  TSOLID(TSOLD,  TLOLD,  KL,  TSNEW) 

IMPLICIT  NONE 
C 

C  VARIABLES  USED  IN  TSOLID 
C 

C  TSOLD  =  VECTOR  OF  POROUS  SOLID  TEMPERATURE  AT  PREVIOUS  TIME  STEP 
C  KL  =  LOCAL  COOLANT  THERMAL  CONDUCTIVITY 
C  DT  =  TIME  STEP 

C  DX  =  X  DISTANCE  STEP  SIZE  THROUGH  WALL 

C  HCL  =  PACKED  BED  OF  SPHERES  HEAT  TRANSFER  COEFFICIENT 

C  CO  =  PROPORTIONALITY  CONSTANT  USED  TO  DETERMINE  SOLID  SURFACE  AREA 

C  R  =  POROUS  SPHERE  RADIUS 

C  EPS  =  POROSITY  OF  THE  POROUS  SOLE) 

C  TLOLD  =  VECTOR  OF  TRANSPIRATION  COOLANT  TEMPERATURE  AT  PREVIOUS 
C  TIME  STEP 

C  TB  =  BULK  TEMPERATURE  OF  THE  COOLANT  IN  THE  COOLANT  CHANNEL 
C  MDOTC  =  MASS  FLOW  RATE  OF  COOLANT  IN  COOLANT  CHANNEL 
C  DH  =  HYDRAULIC  DIAMETER  OF  COOLANT  CHANNEL 
C  PB  =  COOLANT  LINE  PRESSURE  (BULK  CONDITIONS) 

C  PI  =  3.14159 

C  THROATD  =  NOZZLE  AERODYNAMIC  THROAT  DIAMETER 
C  MDOTF  =  MASS  FLOW  RATE  OF  COMBUSTION  PRODUCTS 
C  W  =  VISCOSITY,  TEMPERATURE  EXPONENT 
C  PH  =  HYDRAULIC  PERIMETER  OF  COOLANT  CHANNEL 
C  PBB  =  COOLANT  CHANNEL  PRESSURE  FOR  REGENERATIVE  COOLING  CASE 
C  GAM  =  RATIO  OF  SPECIFIC  HEATS  OF  COMBUSTION  PRODUCTS 
C  TOG  =  STAGNATION  TEMPERATURE  OF  COMBUSTION  PRODUCTS 
C  MDOTL  =  MASS  FLOW  RATE  OF  COOLANT  THROUGH  POROUS  WALL 
C  TG  =  STATIC  TEMPERATURE  OF  COMBUSTION  PRODUCTS  AT  THROAT 
C  TSNEW  =  SOLID  TEMPERATURE  AT  NEW  TIME  STEP 
C  RHOS  =  SOLID  DENSITY 
C  CPS  =  SOLID  SPECIFIC  HEAT 
C  KW  =  SOLID  THERMAL  CONDUCTIVITY 

C  KS  =  EFFECTIVE  THERMAL  CONDUCTIVITY  OF  THE  POROUS  SOLID 
C  HCS  =  COLD  GAS  SIDE  HEAT  TRANSFER  COEFFICIENT 
C  TWH  =  HOT  GAS  SIDE  SOLID  WALL  TEMPERATURE 
C  TAW  =  ADIABATIC  WALL  TEMPERATURE 
C  HG  =  HOT  GAS  SIDE  HEAT  TRANSFER  COEFFICIENT 
C  STABl  =  STABILITY  CONSTRAINT 
C  STAB2  =  STABILITY  CONSTRAINT 
C  STAB3  =  STABILITY  CONSTRAINT 
C  MUOF  =  COMBUSTION  PRODUCTS  VISCOSITY  AT  TOG 
C  CPOF  =  COMBUSTION  PRODUCTS  SPECIFIC  HEAT  AT  TOG 
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C  PROF  =  COMBUSTION  PRODUCTS  PRANDTL  NUMBER  AT  TOG 
C  N  =  INDEX  USED  FOR  SOLID  TEMPERATURE  ARRAY 
C  M  =  INDEX  USED  FOR  COOLANT  TEMPERATURE  ARRAY 
C  LASTX  =  GRID  SIZE 
C 

C  DECLARE  VARIABLES 
C 

DOUBLE  PRECISION  TSOLD,  KL,  DT.  DX,  HCL,  CO,  R,  EPS,  TLOLD,  TB 
DOUBLE  PRECISION  MDOTC,  DH,  PB,  PI,  THROATD,  MDOTF,  W,  PH,  PBB 
DOUBLE  PRECISION  GAM,  TOG,  MDOTL,  TG,  TSNEW,  RHOS,  CPS,  KW 
DOUBLE  PRECISION  KS,  HCS,  TWH,  TAW,  HG,  STABl,  STAB2,  STAB3 
DOUBLE  PRECISION  MUOF,  CPOF,  PROF 

INTEGER  N,  M,  LASTX 

DIMENSION  VARIABLES  THAT  ARE  ARRAYS 

DIMENSION  TSOLD(20000),  KL(20000),  TLOLD(20000) 

DIMENSION  TSNEW(20000) 

THIS  COMMON  BLOCK  CONTAINS  THE  GRID  DIMENSIONS  USED 
COMMON  /DIMEN/  LASTX,  DX 

THIS  COMMON  BLOCK  CONTAINS  OUTPUT  VARIABLES  FROM  TSOLID 
COMMON  /SOLBLK/  HCS,  HG,  TAW 

THIS  COMMON  BLOCK  CONTAINS  VARIABLES  COMMON  TO  TSOLID  AND  TLIQUID 

COMMON  /TSTL/  CO,  EPS,  R,  PB,  TB,  KW,  HCL,  MDOTL 

THIS  COMMON  BLOCK  CONTAINS  VARIABLES  COMMON  TO  TSOLID  ONLY 

COMMON  /TSBLK/  RHOS,  CPS,  MDOTC,  DH,  PI,  THROATD,  MDOTF 
COMMON  /TSBLK/  PH,  W,  GAM,  TOG,  TG,  PBB,  MUOF,  CPOF,  PROF 

BEGIN  THE  MAIN  LOOP  TO  MARCH  THROUGH  THE  POROUS  WALL 

DO  N=1,(LASTX+1) 

SYNCHRONIZE  THE  LIQUID  AND  SOLID  INDICES 
M  =  N  +  2 

CALCULATE  THE  EFFECTIVE  THERMAL  CONDUCTIVITY  OF  THE  SOLID 

KS  =  KL(N)*((2.0*KL(N)+KW)-2.0*(L0-EPS)*(KL(N)-KW))/ 

+  (2.0*KL(N)+KW+(L0-EPS)*(KL(N)-KW)) 

C 
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INITIALIZE  THE  STABILITY  CONSTRAINT  CONSTANTS 


STABI  =  1.0 
STAB2=  1.0 
STAB3  =  1.0 


IF  CALCULATING  THE  COLD  GAS  SIDE  BOUNDARY  TEMPERATURE  USE  THIS 
DF-THEN  LOOP 

IF(N.EQ.1)THEN 

GENERATE  THE  COLD  GAS  SIDE  HEAT  TRANSFER  COEFFICIENT 

CALL  HCOOLSIDE(PBB,  TB,  TSOLD(N),  MDOTC,  DH,  PH,  HCS) 

GENERATE  THE  STABLE  TIME  STEP  FOR  THIS  LOCATION 

DT  =  0.9/((2.0*KS)/((DX**2.0)*RHOS*CPS)+ 

+  (HCL*CO)/(RHOS*R*CPS)+(2.0*HCS)/(RHOS*DX*CPS)) 

CALCULATE  THE  NEW  TEMPERATURE  AT  THE  COLD  GAS  SIDE  BOUNDARY 

TSNEW(N)  =  ((2.0*KS*DT)/((DX**2.0)*RHOS*CPS))*(TSOLD(N+1) 

+  -TSOLD(N))-((HCL*CO*DT)/(RHOS*R*CPS))*(TSOLD(N)- 
+  TLOLD(M))-((2.0*HCS*DT)/(RHOS*DX*CPS))*(TSOLD(N)-TB)+TSOLD(N) 

CALLCULATE  STABILITY  CONSTRAINT 

STABI  =  1.0/((2.0*KS)/((DX**2.0)*RHOS*CPS)+ 

+  (HCL*CO)/(RHOS*R*CPS)+(2.0*HCS)/(RHOS*DX*CPS)) 

END  IF 

IF  CALCULATING  AN  INTERIOR  NODE  POINT  USE  THIS  IF-THEN  LOOP 

IF  ((N.NE.l).AND.(N.NE.(LASTX+l)))  THEN 

CALCULATE  THE  TIME  STEP  FOR  THIS  LOCATION 

DT  =  0.9/((2.0*KS)/((DX**2.0)*RHOS*CPS)+ 

+  (HCL*CO)/(RHOS*R*CPS)) 

CALCULATE  THE  NEW  TEMPERATURE  AT  THIS  LOCATION 

TSNEW(N)  =  ((KS*DT)/((DX**2.0)*RHOS*CPS))*(TSOLD(N+1)+ 

+  TSOLD(N-1)-2.0*TSOLD(N))-((HCL*CO*DT)/(RHOS*R*CPS))* 

+  (TSOLD(N)-TLOLD(M))+TSOLD(N) 

CALCULATE  STABILITY  CONSTRAINT 

STAB2  =  1.0/((2.0*KS)/((DX**2.0)*RHOS*CPS)+ 


63 


+  (HCL*CO)/(RHOS*R*CPS)) 

END  IF 
C 

C  IF  CALCULATING  THE  HOT  GAS  SIDE  BOUNDARY  TEMPERATURE  THEN  USE 
C  THIS  IF-THEN  LOOP 
C 

IF  (N.EQ.(LASTX+I))  THEN 
C 

C  ASSIGN  THE  PREVIOUS  TIME  STEP  TEMPERATURE  TO  THE  HOT  GAS  SIDE 
C  WALL  TEMPERATURE 
C 

TWH  =  TSOLD(N) 

C 

C  GENERATE  THE  HOT  GAS  SIDE  HEAT  TRANSFER  COEFFICIENT 
C 

CALL  HOTGAS(TWH,  THROATD,  MDOTF,  W,  GAM,  TOG,  TG, 

+  MDOTL,  PI,  EPS,  TAW,  MUOF,  CPOF,  PROF,  HG) 

C 

C  CALCULATE  A  TIME  STEP  FOR  THIS  LOCATION 
C 

DT  =  0.9/((2.0*KS)/((DX**2.0)*RHOS*CPS)+ 

+  (HCL*CO)/(RHOS*R*CPS)+(2.0*HG)/(RHOS*DX*CPS)) 

C 

C  CALCULATE  THE  NEW  TEMPERATURE  AT  THIS  LOCATION 
C 

TSNEWCN)  =  ((2.0*HG*DT)/(RHOS*DX*CPS))*(TAW-TSOLD(N))- 
+  ((HCL*CO*DT)/(RHOS*R*CPS))*(TSOLD(N)-TLOLD(M))- 
+  ((2.0*KS*DT)/((DX**2.0)*RHOS*CPS))*(TSOLD(N)-TSOLD(N-I))+TSOLD(N) 

C 

C  CALCULATE  STABILITY  CONSTRAINT 
C 

STAB3  =  1.0/((2.0*KS)/((DX**2.0)*RHOS*CPS)+ 

+  (HCL*CO)/(RHOS*R*CPS)+(2.0*HG)/(RHOS*DX*CPS)) 

END  IF 
C 

C  CHECK  TIME  STEPS  AGAINST  STABILITY  CONSTRAINTS  TO  PREVENT  INSTABILITY 
C  IF  IT  FAILS,  THE  PROGRAM  WRITES  THE  FOLLOWING  OUTPUT  AND  STOPS  THE 
PROGRAM 
C 

IF  ((DT.GT.STABI).OR.(DT.GT.STAB2).OR.(DT.GT.STAB3))  THEN 
WRITE  (*,*)  'STABILITY  CONSTRAINT  VIOLATION' 

WRITE  (*,*)  'STABl  =  ',  STAB  I 
WRITE  (*,*)  'STAB2  =  ',  STAB2 
WRITE  (*,*)  'STAB3  =  ',  STAB3 
WRITE  (*,*)  'DT  USED  =  ',DT 
CLOSE!  14) 

STOP 
END  IF 
C 

C  END  THE  MAIN  LOOP 
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END  DO 
RETURN 

END  TSOLID  SUBROUTINE 
END 

TfflS  SUBROUTINE  CALCULATES  THE  COLD  GAS  SIDE  HEAT  TRANSFER 
COEFFICIENT  GIVEN  COOLANT  BULK  PRESSURE  AND  TEMPERATURE, 

POROUS  SOLID  TEMPERATURE,  MASS  FLOW  RATE  THROUGH  THE  COOLANT 
CHANNEL,  AND  HYDRAULIC  DIAMETER  AND  PERIMETER 

SUBROUTINE  HCOOLSIDE(PB,  TB,  TS,  MDOT,  DH,  PH,  HCS) 

IMPLICIT  NONE 

VARIABLES  USED  IN  HCOOLSIDE 

PB  =  COOLANT  LINE  PRESSURE  (BULK  CONDITIONS) 

TB  =  COOLANT  BULK  TEMPERATURE  IN  COOLANT  CHANNEL 

TS  =  SURFACE  TEMPERATURE  OF  POROUS  SOLID  ON  COLD  GAS  SIDE 

MDOT  =  MASS  FLOW  RATE  OF  COOLANT  IN  COOLANT  CHANNEL 

DH  =  HYDRAULIC  DIAMETER  OF  COOLANT  CHANNEL 

PH  =  HYDRAULIC  PERIMETER  OF  COOLANT  CHANNEL 

HCS  =  COLD  GAS  SIDE  HEAT  TRANSFER  COEFFICIENT 

MUB  =  COOLANT  VISCOSITY  AT  BULK  CONDITIONS 

KB  =  COOLANT  THERMAL  CONDUCTIVITY  AT  BULK  CONDITIONS 

RHOB  =  COOLANT  DENSITY  AT  BULK  CONDITIONS 

CPB  =  COOLANT  SPECIFIC  HEAT  AT  BULK  CONDITIONS 

MULS  =  COOLANT  VISCOSITY  AT  SOLID  SURFACE  TEMPERATURE 

KUS  =  COOLANT  THERMAL  CONDUCTIVITY  AT  SOLID  SURFACE  TEMPERATURE 

RHOLS  =  COOLANT  DENSITY  AT  SOLID  SURFACE  TEMPERATURE 

CPUS  =  COOLANT  SPECIFIC  HEAT  AT  SOLID  SURFACE  TEMPERATURE 

PRB  =  PRANDTL  NUMBER  AT  BULK  CONDITIONS 

REDB  =  REYNOLD'S  NUMBER  AT  BULK  CONDITIONS 

DECLARE  VARIABLES 

DOUBLE  PRECISION  PB,  TB,  TS,  MDOT,  DH,  PH,  HCS 

DOUBLE  PRECISION  MUB,  KB,  RHOB,  CPB,  MULS,  KLS,  RHOLS,  CPLS 

DOUBLE  PRECISION  PRB,  REDB 

GENERATE  HYDROGEN  PROPERTIES  AT  BULK  CONDITIONS 
CALL  H2PROP(PB,  TB,  MUB,  KB,  RHOB,  CPB) 

CALCULATE  PRANDTL  NUMBER  AND  REYNOLD'S  NUMBER 
PRB  =  (MUB*CPB)/KB 
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REDB  =  (4.0*MDOT)/(MUB*PH) 

GENERATE  HYDROGEN  PROPERTIES  AT  SOLID  SURFACE  TEMPERATURE 

CALL  H2PROP(PB,  TS,  MULS,  KLS,  RHOLS.  CPLS) 

CALCULATE  THE  COLD  GAS  SIDE  HEAT  TRANSFER  COEFFICIENT 

HCS  =  ((KB*0.027)/DH)*(REDB**0.8)*(PRB**(LO/3.0))* 

+  ((MUB/MULS)**.I4) 

RETURN 

END  HCOOLSIDE  SUBROUTINE 
END 

THIS  SUBROUTINE  CALCULATES  THE  HOT  GAS  SIDE  HEAT  TRANSFER 
COEFFICIENT  INCLUDING  THE  EFFECT  OF  BLOWING  RATIO 

SUBROUTINE  HOTGAS(TWH,  THROATD,  MOOT,  W,  GAM,  TOG.  TG, 

+  MDOTL,  PL  EPS,  TAW,  MUOF,  CPOF,  PROF,  HG) 

IMPLICIT  NONE 
C 

C  VARIABLES  USED  IN  HOTGAS 
C 

C  TWH  =  HOT  GAS  SIDE  SOLID  SURFACE  TEMPERATURE 
C  THROATD  =  NOZZLE  AERODYNAMIC  THROAT  DIAMETER 
C  MDOT  =  COMBUSTION  GAS  MASS  FLOW  RATE 
C  W  =  VISCOSITY,  TEMPERATURE  EXPONENT 
C  GAM  =  RATIO  OF  SPECIFIC  HEATS  FOR  COMBUSTION  GASES 
C  TOG  =  COMBUSTION  GAS  STAGNATION  TEMPERATURE 
C  TG  =  COMBUSTION  GAS  STATIC  TEMPERATURE  AT  THROAT 
C  MDOTL  =  COOLANT  MASS  FLOW  RATE  THROUGH  POROUS  WALL 
C  PI  =  3.14159 
C  EPS  =  POROSITY 

C  TAW  =  ADIABATIC  WALL  TEMPERATURE 
C  HG  =  HOT  GAS  SIDE  HEAT  TRANSFER  COEFFICIENT 
C  AT  =  NOZZLE  AERODYNAMIC  THROAT  AREA 

C  MUOF  =  COMBUSTION  GAS  VISCOSITY  AT  STAGNATION  TEMPERATURE 
C  CPOF  =  COMBUSTION  GAS  SPECIFIC  HEAT  AT  STAGNATION  TEMPERATURE 
C  BR  =  BLOWING  RATIO 
C  SIGMA  =  BARTZ'  EQUATION  CONSTANT 

C  HGO  =  HOT  GAS  SIDE  HEAT  TRANSFER  COEFFICIENT  WITHOUT  BLOWING 
C  PROF  =  COMBUSTION  GAS  PRANDTL  NUMBER  AT  STAGNATION  TEMPERATURE 
C  RHOFUF  =  MASS  FLUX  OF  COMBUSTION  GASES 
C 

C  DECLARE  VARIABLES 
C 

DOUBLE  PRECISION  TWH,  THROATD,  MDOT,  W,  GAM,  TOG,  TG 
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DOUBLE  PRECISION  MDOTL.  PI.  EPS,  TAW.  HG,  AT,  MUOF,  CPOF 
DOUBLE  PRECISION  BR,  SIGMA,  HGO,  PROF,  RHOFUF 

THIS  COMMON  BLOCK  CONTAINS  OUTPUT  VARIABLES  FROM  HOTGAS 

COMMON  /HGBLK/  BR 

CALCULATE  THE  AERODYNAMIC  THROAT  AREA 
AT  =  (PI*(THROATD**2.0))/4,0 
CALCULATE  THE  ADIABATIC  WALL  TEMPERATURE 
CALL  FINDTAW(TG,  TOG,  TWH,  TAW) 

CALCULATE  THE  MASS  FLUX  OF  THE  COMBUSTION  GASES 
RHOFUF  =  MDOT/AT 
CALCULATE  THE  BLOWING  RATIO 
BR  =  MDOTL/(RHOFUF*EPS) 

THIS  LIMITER  PREVENTS  THE  BLOWING  RATIO  TO  EXCEED  0.01 

IF(BR.GE.1D-02)THEN 
WRITE  (*,*)  'MDOTL  MDOTL 
WRITE  (*,*)  'BR  '.BR 

WRITE  (*,*)  'THE  BLOWING  RATIO  HAS  EXCEEDED  1  PERCENT' 

STOP 
END  IF 

CALCULATE  SIGMA 

SIGMA  =  LO/(((0.5*(TWH/TOG)*((GAM+1.0)/2.0)-H).5)*'»(0.8-0,2*W))* 

+  (((GAM+L0)/2.0)**(0.2*W))) 

CALCULATE  THE  HEAT  TRANSFER  COEFFICIENT  WITHOUT  BLOWING  EFFECTS 

HGO  =  ((0.026/(THROATD**0.2))*(((MU0F**0.2)*CP0F)/(PR0F*'*0.6))* 

+  ((RHOFUF)**0.8)*((2.0)**0.1))*SIGMA 

INCLUDE  BLOWING  EFFECTS 

HG  =  HGO*(L0-38.0*BR) 

RETURN 

END  HOTGAS  SUBROUTINE 
END 
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TfflS  SUBROUTINE  CALCULATES  THE  ADIABATIC  WALL  TEMPERATURE 

SUBROUTINE  FINDTAW(TG,  TOG,  TWH,  TAW) 

IMPLICIT  NONE 

VARIABLES  USED  IN  FINDTAW 

TG  =  STATIC  COMBUSTION  GAS  TEMPERATURE  AT  THROAT 

TOG  =  COMBUSTION  GAS  STAGNATION  TEMPERATURE 

TWH  =  HOT  GAS  SIDE  SOLID  SURFACE  TEMPERATURE 

TAW  =  ADIABATIC  WALL  TEMPERATURE 

DELTA  =  CONVERGENCE  CRITERIA 

TFNEW  =  NEW  FILM  TEMPERATURE 

TFOLD  =  OLD  FILM  TEMPERATURE 

MUF  =  COMBUSTION  GAS  VISCOSITY  AT  FILM  TEMPERATURE 

CPF  =  COMBUSTION  GAS  SPECIFIC  HEAT  AT  FILM  TEMPERATURE 

KF  =  COMBUSTION  GAS  THERMAL  CONDUCTIVITY  AT  FILM  TEMPERATURE 

PRF  =  COMBUSTION  GAS  PRANDTL  NUMBER  AT  FILM  TEMPERATURE 

DECLARE  VARIABLES 

DOUBLE  PRECISION  TG,  TOG,  TWH,  TAW,  DELTA,  TFNEW,  TFOLD,  MUF 
DOUBLE  PRECISION  CPF,  KF,  PRF 

INITIALIZE  VARIABLES 

DELTA  =  9999,0 
TFOLD  =  (TWH+T0G)/2.0 

BEGIN  MAIN  LOOP  AND  ITERATE  UNTIL  CONVERGED  ON  CORRECT  FILM 
TEMPERATURE 

DO  WHILE  (DELTA.GE.lD-04) 

GENERATE  COMBUSTION  GAS  PROPERTIES  AT  OLD  FILM  TEMPERATURE 

CALL  HOTGASPROP(TFOLD,  MUF,  CPF,  KF) 

CALCULATE  NEW  PRANDTL  NUMBER,  ADIABATIC  WALL  TEMPERATURE  AND 
FILM  TEMPERATURE 

PRF  =  (MUF*CPF)/KF 

TAW  =  (PRF**(L0/3.0))*(T0G-TG)+TG 

TFNEW  =  TWH-K).23*(TG-TWH)-H).I9*(TAW-TWH) 

COMPARE  OLD  AND  NEW  FILM  TEMPERATURES  TO  CHECK  FOR  CONVERGENCE 

DELTA  =  (TFNEW-TFOLD)**2.0 
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TFOLD  =  TFNEW 


END  MAIN  LOOP 

END  DO 
RETURN 

END  FINDTAW  SUBROUTINE 
END 

TfflS  SUBROUTINE  GENERATES  THE  COMBUSTION  GAS  PROPERTIES  AT 
DIFFERENT  TEMPERATURES 

SUBROUTINE  HOTGASPROP(T,  MUF,  CPF,  KF) 

IMPLICIT  NONE 

VARIABLES  USED  IN  HOTGASPROP 

T  =  LOCAL  STATIC  TEMPERATURE  AT  THROAT 

MUF  =  COMBUSTION  GAS  VISCOSITY 

CPF  =  COMBUSTION  GAS  SPECIFIC  HEAT 

KF  =  COMBUSTION  GAS  THERMAL  CONDUCTIVITY 

PRF  =  COMBUSTION  GAS  PRANDTL  NUMBER 

DECLARE  VARIABLES 

DOUBLE  PRECISION  T,  MUF,  CPF,  KF,  PRF 

GENERATE  THERMODYNAMIC  PROPERTIES 

MUF  =  (-3.19123D-12)*(T**2.0)+(4.21474D-08)*T-(L9I043D-06) 

CPF  =  (-1.7605 1D-13)*(T**5.0)+(2.389596D-09)*(T**4.0) 

+  -(1. 12842D-05)*(T**3.0)+(0.0240347)*(T**2.0)-(22.7384*T)+ 

+  10438 

PRF  =  (-3.63422D-27)*(T**8,0)+(6,93638D-23)*(T**7.0)- 
+  (5.70795D-19)*(T**6.0)+(2.63999D-15)*(T**5.0)- 
+  (7,47121D-12)*(T**4.0)+(L31513D-8)*(T**3.0)- 
+  (1,3  9728D-5)*(T**2.0)+(0.008 1 8492)*T-L4202 
KF  =  MUF*CPF/PRF 
RETURN 

END  HOTGASPROP  SUBROUTINE 
END 

THIS  SUBROUTINE  SENDS  THE  TEMPERATURE  PROFILES  TO  THE  SPECIFIED 
FILENAME,  AND  THE  SCREEN 
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SUBROUTINE  OUTPUT(TSNEW,  TLNEW) 


IMPLICIT  NONE 
C 

C  VARIABLES  USED  IN  OUTPUT 
C 

C  DX  =  INCREMENTAL  DISTANCE  STEP 

C  TSNEW  =  SOLID  TEMPERATURE  PROFILE  AT  NEW  TIME  STEP 
C  TLNEW  =  COOLANT  TEMPERATURE  PROFILE  AT  NEW  TIME  STEP 
C  PB  =  COOLANT  LINE  PRESSURE 
C  R  =  POROUS  SPHERE  RADIUS 

C  CO  =  PROPORTIONALITY  CONSTANT  FOR  POROUS  SOLID  SURFACE  AREA 
C  EPS  =  POROSITY 

C  KW  =  SOLID  THERMAL  CONDUCTIVITY 

C  MDOTL  =  COOLANT  MASS  FLOW  RATE  THROUGH  POROUS  SOLID 
C  BR  =  BLOWING  RATIO 

C  HCS  =  COLD  GAS  SIDE  HEAT  TRANSFER  COEFFICIENT 
C  HG  =  HOT  GAS  SIDE  HEAT  TRANSFER  COEFFICIENT 
C  HCL  =  PACKED  BED  OF  SPHERES  HEAT  TRANSFER  COEFFICIENT 
C  TAW  =  ADIABATIC  WALL  TEMPERATURE 
C  TB  =  COOLANT  BULK  TEMPERATURE 
C  N  =  INDEX  USED  FOR  SOLID  TEMPERATURE  ARRAY 
C  M  =  INDEX  USED  FOR  COOLANT  TEMPERATURE  ARRAY 
C  LASTX  =  GRID  SIZE 
C  FILE  1  =  OUTPUT  FILENAME 
C 

C  DECLARE  VARIABLES 
C 

DOUBLE  PRECISION  DX,  TSNEW,  TLNEW,  PB,  R,  CO,  EPS,  KW 
DOUBLE  PRECISION  MDOTL,  BR,  HCS,  HG,  HCL,  TAW,  TB 

INTEGER  N,  M,  LASTX 

CHARACTER*64  FILE  I 

DIMENSION  VARIABLES  THAT  ARE  ARRAYS 
DIMENSION  TSNEW(20000),  TLNEW(20000) 

THIS  COMMON  BLOCK  CONTAINS  THE  GRID  DIMENSIONS 
COMMON  /DIMEN/  LASTX,  DX 

THIS  COMMON  BLOCK  CONTAINS  VARIABLES  COMMON  TO  TSOLID  AND  TLIQUID 
COMMON  /TSTL/  CO,  EPS,  R,  PB,  TB,  KW,  HCL,  MDOTL 
THIS  COMMON  BLOCK  CONTAINS  THE  OUTPUT  FILENAME 
COMMON  /FILEBLK/  FILE  I 
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TfflS  COMMON  BLOCK  CONTAINS  OUTPUT  DATA  FROM  HOTGAS 
COMMON  /HGBLK/  BR 

THIS  COMMON  BLOCK  CONTAINS  OUTPUT  DATA  FROM  TSOLID 

COMMON  /SOLBLK/  HCS,  HG,  TAW 

OPEN  THE  OUTPUT  FE.E 

OPEN  (UNIT=14,  FILE=FILE1) 

WRITE  THE  DATA  TO  THE  OUTPUT  FILE 

DO  N=L(LASTX+1) 

M  =  N+2 

WRITE  (14,10)  1000.0*(N-1)*DX,  TSNEW(N),  TLNEW(M) 

END  DO 

CLOSE(I4) 

WRITE  THE  INPUT  DATA  AND  OTHER  VARIABLES  TO  THE  SCREEN 

WRITE  (*.*)  'INPUT  VARIABLES' 

WRITE  (*,*) 

WRITE  (*,*)  'PB  =  ',PB 
WRITE  (*,*)  'R  =  ',R 
WRITE  (*,*)  'CO  =  ’,C0 
WRITE  (*,*) 'EPS  =',EPS 
WRITE  (*,*)  'KW  =  ',KW 
WRITE  (*,*) 

WRITE  (*,*)  'DATA  IS  STORED  IN  FILE  =  ',FILE1 
WRITE  (*,*) 

WRITE  (*,*)  'MDOTL  =  ',MDOTL 
WRITE  (*,♦)  'BR  =  ',BR 
WRITE  (*,*) 'HCS  =',HCS 
WRITE  (*,*) 'HCL  =',HCL 
WRITE  (*,*)  'HG  =  ',HG 
WRITE  (*,*) 'TAW  =',TAW 

FORMATTED  OUTPUT  STATEMENT 

)  FORMAT  (IX,FI0.8,IX,F15.9,1X,F15.9) 

RETURN 

END  OUTPUT  SUBROUTINE 

END 
C 
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TfflS  SUBROUTINE  CALCULATES  THE  COOLANT  TEMPERATURE  PROFILE 

SUBROUTINE  TLIQUID(TSOLD,  TLOLD,  KL,  MULAVG,  RHOLAVG,  TLNEW) 

IMPLICIT  NONE 
C 

C  VARIABLES  USED  IN  TLIQUID 
C 

C  TLOLD  =  COOLANT  TEMPERATURE  PROFILE  AT  PREVIOUS  TIME  STEP 
C  MDOTL  =  COOLANT  MASS  FLOW  RATE  THROUGH  POROUS  SOLID 
C  EPS  =  POROSITY 
C  R  =  POROUS  SPHERE  RADIUS 
C  TB  =  COOLANT  BULK  TEMPERATURE 
C  DX  =  INCREMENTAL  DISTANCE  STEP  SIZE 
C  CO  =  PROPORTIONALITY  CONSTANT  FOR  SOLID  SURFACE  AREA 
C  TSOLD  =  SOLID  TEMPERATURE  PROFILE  AT  PREVIOUS  TIME  STEP 
C  KL  =  COOLANT  THERMAL  CONDUCTIVITIES 
C  RHOL  =  COOLANT  DENSITIES 
C  MULAVG  =  AVERAGE  COOLANT  VISCOSITY 
C  RHOLAVG  =  AVERAGE  COOLANT  DENSITY 
C  TLNEW  =  COOLANT  TEMPERATURE  PROFILE  AT  NEW  TIME  STEP 
C  KW  =  SOLID  THERMAL  CONDUCTIVITY 
C  DELTAP  =  PRESSURE  DIFFERENCE  ACROSS  POROUS  WALL 
C  PERMEABLE  =  PERMEABILITY 
C  PL  =  LOCAL  COOLANT  PRESSURE 
C  KLL  =  LOCAL  COOLANT  THERMAL  CONDUCTIVITY 
C  CPL  =  COOLANT  SPECIFIC  HEATS 
C  UL  =  COOLANT  VELOCITIES 

C  HCL  =  PACICED  BED  OF  SPHERES  HEAT  TRANSFER  COEFFICIENT 
C  PB  =  COOLANT  CHANNEL  PRESSURE 
C  DP  =  INCREMENTAL  PRESSURE  DIFFERENCE 
C  Cl,  C2,  C3  =  PLACE  HOLDER  CONSTANTS 

C  PLAVG  =  AVERAGE  COOLANT  PRESSURE  WITHIN  POROUS  WALL 
C  TLAVG  =  AVERAGE  COOLANT  TEMPERATURE  WITHIN  POROUS  WALL 
C  RHOB  =  COOLANT  DENSITY  AT  BULK  CONDITIONS 
C  THICKNESS  =  POROUS  WALL  THICKNESS 
C  MULL  =  LOCAL  COOLANT  VISCOSITY 
C  CPLL  =  LOCAL  COOLANT  SPECIFIC  HEAT 
C  ULL  =  LOCAL  COOLANT  VELOCITY 
C  CPLAVG  =  AVERAGE  COOLANT  SPECIFIC  HEAT 
C  KLAVG  =  AVERAGE  COOLANT  THERMAL  CONDUCTIVITY 
C  N  =  INDEX  USED  FOR  SOLID  TEMPERATURE  ARRAY 
C  M  =  INDEX  USED  FOR  COOLANT  TEMPERATURE  ARRAY 
C  LASTX  =  GRID  SIZE 
C 

C  DECLARE  VARIABLES 
C 

DOUBLE  PRECISION  TLOLD,  MDOTL,  EPS,  R,  TB,  DX,  CO,  TSOLD,  KL,  RHOL 
DOUBLE  PRECISION  MULAVG,  RHOLAVG,  TLNEW,  KW,  DELTAP,  PERMEABLE 
DOUBLE  PRECISION  PL,  KLL,  RHOLL,  CPL,  UL,  HCL,  PB,  DP 
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DOUBLE  PRECISION  Cl.  C2,  C3,  PLAVG,  TLAVG,  RHOB.  THICKNESS 
DOUBLE  PRECISION  MULL,  CPLL,  ULL,  CPLAVG,  KLAVG 

INTEGER  N,  M,  LASTX 

DIMENSION  VARIABLES  THAT  ARE  ARRAYS 

DIMENSION  TLOLD(20000),  TSOLD(20000),  KL(20000) 

DIMENSION  RHOL(20000),  CPL(20000),  UL(20000) 

DIMENSION  TLNEW(20000) 

THIS  COMMON  BLOCK  CONTAINS  THE  GRID  DIMENSIONS 
COMMON  /DIMEN/  LASTX,  DX 

THIS  COMMON  BLOCK  CONTAINS  VARIABLES  COMMON  TO  TSOLID  AND  TLIQUID 

COMMON  /TSTL/  CO,  EPS,  R,  PB,  TB,  KW,  HCL,  MDOTL 

THIS  COMMON  BLOCK  CONTAINS  VARIABLES  COMMON  TO  TLIQUID  ONLY 

COMMON  /TLBLK/  PLAVG,  TLAVG,  RHOB,DP,DELTAP,PERMEABLE,THICKNESS 

CALCULATE  THE  COOLANT  MASS  FLOW  RATE  THROUGH  THE  POROUS  WALL 

MDOTL  =  (DELTAP*PERMEABLE*RHOLAVG)/(THICKNESS*MULAVG) 

CALCULATE  THE  PACKED  BED  OF  SPHERES  HEAT  TRANSFER  COEFFICIENT 

CALL  HCOOL(PLAVG,  TLAVG,  MDOTL,  RHOB,  EPS,  R,  HCL) 

BEGIN  MAIN  LOOP  TO  STEP  THROUGH  THE  POROUS  WALL 

DON=l,(LASTX+l) 

M=N  +  2 

CALCULATE  THE  LOCAL  COOLANT  PRESSURE 
PL  =  PB-(DP*(N-1)) 

GENERATE  HYDROGEN  PROPERTIES  AT  THE  LOCAL  TEMPERATURE 

CALL  H2PROP(PL,  TLOLD(M),  MULL,  KLL,  RHOLL,  CPLL) 

KL(N)  =  KLL 
RHOL(N)  =  RHOLL 
CPL(N)  =  CPLL 

CALCULATE  THE  LOCAL  COOLANT  VELOCITY 
ULL  =  MDOTL/(EPS*RHOLL) 
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UL(N)  =  ULL 


SET  THE  FIRST  TWO  TEMPERATURES  BEFORE  THE  POROUS  WALL  EQUAL 
TO  THE  BULK  TEMPERATURE 

TLNEW(1)  =  TB 
TLNEW(2)  =  TB 

DEFINE  THE  PLACEHOLDER  CONSTANTS  Cl  AND  C2 

Cl  =  KLL/(DX**2.0) 

C2  =  (RHOLL*ULL*CPLL)/DX 

DEFINE  THE  PLACEHOLDER  CONSTANT  C3  BASED  UPON  FIRST  OR  LAST  GRID 
POINTS 

IF  ((N.NE.  1).AND.(N.NE.(LASTX+1)))  THEN 
C3  =  (HCL+CO)/(EPS*R) 

ELSE 

C3  =  (HCL*C0)/(2.0*EPS*R) 

END  IF 

CALCULATE  THE  NEW  COOLANT  TEMPERATURE  AT  THIS  LOCATION 

TLNEW(M)  =  (1.0/(C2+C3-C1))*(C1*(TLNEW(M.2)- 
+  2.0*TLNEW(M-1))+C2*TLNEW(M-1)  +  C3*TSOLD(N)) 

END  THE  MAIN  LOOP 

END  DO 

CALCULATE  THE  AVERAGE  TEMPERATURE 
TLAVG  =  (TLNEW(3)  +  TLNEW(LASTX+3))/2.0 

GENERATE  HYDROGEN  PROPERTIES  AT  THE  AVERAGE  PRESSURE  AND  TEMPERATURE 

CALL  H2PROP(PLAVG,  TLAVG,  MULAVG,  KLAVG,  RHOLAVG,  CPLAVG) 

RETURN 

END  TLIQUID  SUBROUTINE 
END 

THIS  SUBROUTINE  CALCULATES  THE  PACIOED  BED  OF  SPHERES  HEAT  TRANSFER 
COEFFICIENT 

SUBROUTINE  HCOOL(PLAVG,  TLAVG,  MDOTL,  RHOB,  EPS,  R,  HCL) 

IMPLICIT  NONE 
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C 

C  VARIABLES  USED  IN  HCCXDL 
C 

C  PLAVG  =  AVERAGE  COOLANT  PRESSURE  THROUGH  POROUS  WALL 
C  TLAVG  =  AVERAGE  COOLANT  TEMPERATURE  THROUGH  POROUS  WALL 
C  MDOTL  =  COOLANT  MASS  FLOW  RATE  THROUGH  POROUS  WALL 
C  PR  =  COOLANT  PRANDTL  NUMBER 
C  MU  =  COOLANT  VISCOSITY 
C  K  =  COOLANT  THERMAL  CONDUCTIVITY 
C  RHO  =  COOLANT  DENSITY 
C  CP  =  COOLANT  SPECIFIC  HEAT 
C  EPS  =  POROSITY 
C  R  =  POROUS  SPHERE  RADIUS 

C  HCL  =  PACKED  BED  OF  SPHERES  HEAT  TRANSFER  COEFFICIENT 
C  RHOB  =  COOLANT  DENSITY  AT  BULK  CONDITIONS 
C  REDP  =  COOLANT  REYNOLD'S  NUMBER 
C 

C  DECLARE  VARIABLES 
C 

DOUBLE  PRECISION  PLAVG,  TLAVG,  MDOTL,  PR,  MU,  K,  RHO,  CP 
DOUBLE  PRECISION  EPS,  R,  HCL,  RHOB,  REDP 

GENERATE  HYDROGEN  PROPERTIES  AT  AVERAGE  PRESSURE  AND  TEMPERATURE 
CALL  H2PROP(PLAVG,  TLAVG,  MU,  K,  RHO,  CP) 

CALCULATE  PRANDTL  NUMBER  AND  REYNOLD’S  NUMBER 
PR  =  (MU*CP)/K 

REDP  =  (2,0*R*MDOTL)/(MU*EPS*(1.0-EPS)) 

CALCULATE  THE  PACKED  BED  OF  SPHERES  HEAT  TRANSFER  COEFFICIENT 

HCL  =  (K/(2.0*R))*((1.0-EPS)/EPS)*(0.5*(REDP**0.5)+ 

+  (0.2*(REDP**(2.0/3.0))))*(PR**(LO/3.0)) 

RETURN 

END  HCOOL  SUBROUTINE 
END 
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